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QQ , As demonstrated by Slepian et. al. in a sequence of classical papers (see [32], |33| . [14] . 

[34] . [35]). prolate spheroidal wave functions (PSWFs) provide a natural and efficient tool for 
computing with bandlimited functions defined on an interval. Recently, PSWFs have been 
■^>i ■ becoming increasingly popular in various areas in which such functions occur - this includes 

j /^ physics (e.g. wave phenomena, fluid dynamics), engineering (signal processing, filter design), 

etc. 

To use PSWFs as a computational tool, one needs fast and accurate numerical algorithms for 
the evaluation of PSWFs and related quantities, as well as for the construction of corresponding 
quadrature rules, interpolation formulas, etc. During the last 15 years, substantial progress has 
been made in the design of such algorithms - see, for example, [3Z] (see also [3J, [53], [2], [31] 
for some classical results). 

The complexity of many of the existing algorithms, however, is at least quadratic in the 
band limit c. For example, the evaluation of the nth eigenvalue of the prolate integral operator 
requires at least 0(c 2 ) operations (see e.g. [37]); the construction of accurate quadrature rules 
for the integration (and associated interpolation) of bandlimited functions with band limit c 
requires 0(c 3 ) operations (see e.g. [4]). Therefore, while the existing algorithms are satisfactory 
for moderate values of c (e.g. c < 10 3 ), they tend to be relatively slow when c is large (e.g. 
c>10 4 ). 

In this paper, we describe several numerical algorithms for the evaluation of PSWFs and 
related quantities, and design a class of PSWF-based quadratures for the integration of ban- 
dlimited functions. While the analysis is somewhat involved and will be published separately 
(currently, it can be found in [25], [26]), the resulting numerical algorithms are quite simple and 
k>( ■ efficient in practice. For example, the evaluation of the nth eigenvalue of the prolate integral 

j_j | operator requires 0(n + c-log c) operations; the construction of accurate quadrature rules for the 

integration (and associated interpolation) of bandlimited functions with band limit c requires 
0(c) operations. All algorithms described in this paper produce results essentially to machine 
precision. Our results are illustrated via several numerical experiments. 
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1 Introduction 



The principal purpose of this paper is to describe several numerical algorithms associated with 
bandlimited functions. While these algorithms are quite simple and efficient in practice, the analysis 
is somewhat involved, and will be published separately (currently the proofs and additional details 
can be found in [H], [37], [IB])- 

A function / : K — > R is said to be bandlimited with band limit c > if there exists a function 
a e L 2 [— 1, 1] such that 



In other words, the Fourier transform of a bandlimited function is compactly supported. While ([T]) 
defines / for all real x, one is often interested in bandlimited functions whose argument is confined 
to an interval, e.g. — 1 < x < 1. Such functions are encountered in physics (wave phenomena, fluid 
dynamics), engineering (signal processing), etc. (see e.g. [32], [7J, [2"5]). 

About 50 years ago it was observed that the eigenfunctions of the integral operator F c : L 2 [— 1, 1] — > 
L 2 [—1,1], defined via the formula 



provide a natural tool for dealing with bandlimited functions defined on the interval [—1,1]. More- 
over, it was observed (see [33], [T3], [31]) that the eigenfunctions of F c are precisely the prolate 
spheroidal wave functions (PSWFs), well known from the mathematical physics (see, for example, 



Obviously, to use PSWFs as a computational tool, one needs fast and accurate numerical al- 
gorithms for the evaluation of PSWFs and related quantities, as well as for the construction of 
quadratures, interpolation formulas, etc. For the last 15 years, substantial progress has been made 
in the design of such algorithms - see, for example, [37] (see also [3J, [33], [13], [33] for some classical 
results). 

The complexity of many of the existing algorithms, however, is at least quadratic in the band 
limit c. For example, the evaluation of the nth eigenvalue of the prolate integral operator requires 
0{c 2 + n 2 ) operations (see e.g. [37]); also, the construction of accurate quadrature rules for the 
integration (and associated interpolation) of bandlimited functions with band limit c requires 0(c 3 ) 
operations (see e.g. [3])- Therefore, while the existing algorithms are satisfactory for moderate 
values of c (e.g. c < 10 3 ), they tend to be relatively slow when c is large (e.g. c > 10 4 ). 

In this paper, we describe several numerical algorithms for the evaluation of PSWFs and related 
quantities, and design a class of PSWF-based quadratures for the integration of bandlimited func- 
tions. While the analysis is somewhat involved and will be published separately (currently, it can be 
found in |25) . [26] ). the resulting numerical algorithms are quite simple and efficient in practice. For 
example, the evaluation of the nth eigenvalue of the prolate integral operator requires 0(n + clogc) 
operations; also, the construction of accurate quadrature rules for the integration of bandlimited 
functions with band limit c requires 0{c) operations. In addition, the evaluation of the nth PSWF 
is done in two steps. First, we carry out a certain precomputation, that requires 0(n + clogc) 
operations. Then, each subsequent evaluation of this PSWF at a point in [—1, 1] requires 0(1) 
operations. 

This paper is organized as follows. Section [2J contains a brief overview. Section [3J contains 
mathematical and numerical preliminaries to be used in the rest of the paper. Section [4] contains 
the summary of the principal analytical results of the paper. Section [5] contains the description and 
analysis of the numerical algorithms for the evaluation of the quadrature rules and some related 
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quantities. In Section [SJ we report some numerical results. In Section we illustrate the analysis 
via several numerical experiments. 

2 Overview 

In this section, we provide an overview of the paper. More specifically, Section |2"7T1 is dedicated to the 
numerical evaluation of PSWFs and related quantities. In Section 12.21 we discuss several existing 
quadrature rules for the integration of bandlimited functions. In Section 12.31 we introduce a new 
class of PSWFs-based quadrature rules and describe the underlying ideas. In Section l2~4l we outline 
the analysis (further details can be found in [25 , [26 .). 

2.1 Numerical Evaluation of PSWFs 

For any real c > and integer n > 0, the corresponding PSWF tp n can be expanded into an 
infinite series of Legendre polynomials (see Section |3~2")) . The coefficients of such expansions decay 
superalgebraically (see e.g [37]); in particular, relatively few terms of the Legendre series are required 
to evaluate ip n (x) to essentially the machine precision, for any — 1 < x < 1. The use of this 
observation for the numerical evaluation of PSWFs goes back at least to the classical Bouwkamp 
algorithm [3J (see also Section l3T2l in particular TheoremHUland Remark^ and [37] for more details). 

Thus, the evaluation of PSWFs reduces to the evaluation of the corresponding Legendre coeffi- 
cients. For any integer n > 0, the Legendre coefficients of all the first n PSWFs if)o,ipi, . . . ,ipn-i 
can be obtained via the solution of a certain symmetric tridiagonal eigenproblem roughly of order 
max {n,c} (see Theorem ITUl and Remark [H] in Section I3T1Z1 and also [37 for more details about this 
algorithm). The corresponding eigenvalues Xoj Xi> • • - j Xn-i of the prolate differential operator (see 
Theorem [3] in Section 13- If) are obtained as a by-product of this procedure. On the other hand, 
additional computations are required to evaluate the corresponding eigenvalues Ao, Ai, . . . , \ n -i of 
the integral operator F c (see ^ in Section [TJ. In practice, it is sometimes desirable to evaluate 
extremely small A/s (e.g. 1E-50), which presents a numerical challenge (see Section I3TT|) . To over- 
come this obstacle, the algorithm of [37] evaluates Ao, Ai, . . . , A„_i by computing the ratios A^/Aj+i, 
which turns out to be a well-conditioned numerical procedure (see |37j for more details) . 

Suppose, on the other hand, that one is interested in a single PSWF ip n only (as opposed to all 
the first n PSWFs). Obviously, one can use the algorithm of [37]; however, its cost is at least 0(n 2 ) 
operations (see Remark Moreover, the cost of evaluating the corresponding eigenvalue A„ of the 
prolate integral operator F c (see @) via the algorithm of [37] is at least 0(n 2 ) operations, with a 
large proportionality constant. 

In this paper, we describe more efficient algorithms for the numerical evaluation of ip n and 
associated quantities. In particular, the cost of the evaluation of the Legendre coefficients of ijj n via 
this algorithm is 0(n + clogc) operations (see Section ISTTj) . In addition, the cost of the evaluation 
of the eigenvalue A„ is also 0{n + clogc) operations (see Section I5T2]) . On the other hand, this 
algorithm has the same accuracy as that of [37] : in other words, all of the quantities are evaluated to 
essentially the machine precision (see Section [S] for more details). Since A„ can be extremely small, 
the fact that it can be evaluated to high relative accuracy (without computing the preceding Aj's) 
is, perhaps, surprising (the related analysis is somewhat subtle, and will be published separately; 
see [37J, [35] for some preliminary results). 

2.2 Quadrature Rules for Bandlimited Functions 

One of principal goals of this paper is a class of quadrature rules designed for the integration of 
bandlimited functions with a specified band limit c > over the interval [—1, 1]. Suppose that n > 



3 



is an integer; a quadrature rule of order n is a pair (t^ 1 
of length n, where 



An) w (n) 



-l < 4 n) < 

are referred to as "the quadrature nodes", and 



(n) 



< ti n) < 1 



w {r - 

i ~ ' n 



. Wi n) ) of finite sequences 



(3) 



(4) 



are referred to as "the quadrature weights". If / : [— 1, 1] — > M is a bandlimited function (see ((T|) in 
Section [1}, we use the quadrature rule to approximate the integral of / over the interval [— 1, 1] by 
a finite sum; more specifically, 



L n, 
1 .7 = 1 



(5) 



The PSWFs constitute a natural basis for the bandlimited functions with band limit c > over the 
interval [—1,1] (see Section[TJabove). Therefore, when designing a quadrature rule for the integration 
of such functions, it is reasonable to require that this quadrature rule integrate several first PSWFs 
with band limit c to high accuracy. To describe this property in a more precise manner, we introduce 
the following definition. 

Definition 1. Suppose that c > is a real number, and that n > is an integer. Suppose also that a 
quadrature rule for the integration of bandlimited functions with band limit c over [—1, 1] is specified 
via its n nodes and weights, as in (|5|). Q. Suppose furthermore that e > is a real number, and 
that this quadrature rule integrates the first n PSWFs of band limit c to precision e, in other words, 



^ m (t)dt-j2w^ n) A 

3 = 1 



(*?>) 



(6) 



for every integer m = 0, l,...,n— 1, where tjj rn : [—1,1] -> 1 is the mth PSWF corresponding 
to band limit c. We refer to such quadrature rules as "quadrature rules of order n to precision e 
(corresponding to band limit c)". We omit the reference to c whenever the band limit is clear from 
the context. 

Remark 1. Obviously, if e is the machine precision (e.g. e ~ 1D-16 in double precision calcula- 
tions ), then quadrature rules of order n to precision e (in the sense of Definition [TJ) integrate the 
first n PSWFs exactly, for all practical purposes. 

Remark 2. In practice, for a quadrature rule of order n to precision e to be of any use for the 
integration of bandlimited functions with band limit c, not only e should be "small", but also n has 
to be at least equal to 2c/ir. See Section \3.1\ and f37V for more details. 

Quadrature rules for the integration of bandlimited functions have already been discussed in the 
literature, for example: 

Generalized Gaussian Quadrature Rules. Suppose that n > is an integer, and that 
fx , /2, . . . , f2n are 2n linearly independent functions defined on an interval. Under very mild con- 
ditions on fx, . . . , f2 n , there exists a quadrature rule of order n that integrates these 2n functions 
exactly; moreover, its weights are usually positive. Such quadrature rules are referred to as "gener- 
alized Gaussian quadrature rules" , and their existence was first observed more than 100 years ago 
(see, for example, [12] . |13j . [17] . [TS]). Perhaps surprisingly, numerical algorithms for the design 
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of generalized Gaussian quadrature rules were constructed only recently (see, for example, [3], |16| . 
[55]). These algorithms tend to be rather expensive (they require 0(n 3 ) operations with a large 
proportionality constant). Thus, the evaluation of the nodes and weights of a PSWF-based gener- 
alized Gaussian quadrature rule for accurate integration of bandlimited functions with band limit c 
requires 0(c 3 ) operations (see Remark [2] above, and also [37] for more details). 

Remark 3. We observe that a PSWF-based generalized Gaussian quadrature rule of order n in- 
tegrates the first 2n PSWFs exactly; in other words, (O holds for every integer m between and 
2n — 1 with e = 0. 

Quadrature Rules from [37J. Suppose now that n > is an integer, and that tp n is the nth 
PSWF corresponding to band limit c. Suppose also that t\, . . . , t n are the roots of ip n in the interval 
(— 1, 1) (see Theorem[T]in Section l3~Tj) . Suppose furthermore that Wi, . . . , W n are real numbers, and 
that 

n -i 
53^m(ti)-Wi= / lp m (t)dt, (7) 
z=l 

for every m = 0, . . . , n — 1. Obviously, due to ([7|). the quadrature rule with nodes ti, . . . , t n and 
weights Wi, . . . , W n integrates the first n PSWFs exactly (i.e. ^ holds for every m = 0, . . . , n — 1 
with e = 0). While this quadrature rule is clearly "sub-optimal" compared to the generalized 
Gaussian quadrature rule of order n (the latter integrates the first 2n PSWFs exactly) , it is somewhat 
less expensive to evaluate. More specifically, the cost of evaluating the roots ti, . . . , t n of ip n in (—1, 1) 
and the weights W\, . . . , W n , defined via (J7J|, is dominated by the cost of solving the dense n by 
n linear system (J7J for the unknowns W\ , . . . , W n (see [37] for more details about the numerical 
aspects of this procedure). Thus, due to Remark [5] above, the cost of evaluating the nodes and 
weights of this quadrature rule for accurate integration of bandlimited functions with band limit c 
requires 0(c 3 ) operations. 

Remark 4. The cost of the evaluation of the quadrature rule, defined via {7J), is 0(c 3 ) operations. 
The cost of the evaluation of the generalized Gaussian quadrature rule is also 0(c 3 ) operations, but 
tends to have a larger proportionality constant. 

Remark 5. The quadrature rule defined via (J7J| is based on the PSWFs corresponding to band limit 
c. It turns out, however, that this quadrature rule will also integrate bandlimited functions with band 
limit 2c to high accuracy. The reason for this is that the classical Euclid algorithm for polynomial 
division can be generalized to the PSWFs; the reader is referred to JffTf for further details. 

In this paper, we describe another class of quadrature rules whose nodes are the n roots of ip n 
in ( — 1, 1). However, their weights differ slightly from those defined via ([7]). In particular, strictly 
speaking, these quadrature rules do not integrate the first n PSWFs exactly, as opposed to the 
generalized Gaussian quadrature rules and those defined via ([7]) above. Nevertheless, for any e > 0, 
they do integrate the first n PSWFs to precision e, provided that 

n> — + 10 + A. (log c ). log! (8) 
7T ir z £ 

(see Theorem [j~5l from Section 1431 and Conjectures [3j [4] from Section [7] for more precise statements, 
and Experiment 3 in Section [7. II for some numerical results). 

Thus, provided that e is the machine precision and that (jSJ holds, the quadrature rules of this 
paper are, for all practical purposes, as accurate as those defined via ([TJ above. Also, their nodes 
and weights can be used as starting points for an iterative scheme that computes the generalized 
Gaussian quadrature rule (see, for example, [4], [16], [38] for more details). Last but not least, 
the quadrature rules of this paper are much faster to evaluate than those described above: 0{c) 
operations are required (see Sections 15.31 15.41) . 
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2.3 Intuition Behind Quadrature Weights 



In this section, we describe the quadrature rules of this paper, and discuss the intuition behind them. 

We start with a classical interpolation problem. Suppose that t\ , . . . , t n are n distinct points on 
the interval (—1,1). We need to find the real numbers W\, . . . , W n such that 



/l n 
- 1 i=i 



(9) 



for all polynomials p of degree at most n — 1. In other words, the quadrature rule with nodes 
t\ , . . . , t n and weights W\ , ■ ■ ■ , W n integrates all polynomials of degree up to n — 1 exactly (see ([3]) , 

©, ©). 

To this end, one constructs n polynomials li, . . . , l n of degree n — 1 with the property 

[1 1 = 3 

for every integer i,j = l,...,n (see, for example, |11|). It is easy to verify that, for every j = 1, . . . , n, 
the polynomial lj is defined via the formula 



l.( t ) = ^ 



(11) 



3 J \" "J J 

for all real — 1 < t < 1, where u>„ is defined via the formula 



w n (t) = (t-h) ■ (t-t 2 ) (t-tn), (12) 

for all real — 1 < t < 1 (in other words, w n is the polynomial of degree n whose roots are precisely 
ti,...,t n ). The weights W\, . . . , W„ are defined via the formula 

for every integer j = 1, . . . , n. 

In our case, the basis functions are the PSWFs rather than polynomials. We will consider 
the quadrature rule (t%, . . . ,t n , Wi, . . . , W n ), with t\, . . . , t n the roots of i\) n on the interval (—1,1), 
and Wi, ■ ■ ■ , W n to be determined. If we choose the weights W\,..., W n such that the resulting 
quadrature rule integrates the first n PSWFs exactly, this will lead to the linear system ([7]) from 
Section [2.21 (and hence to the corresponding quadrature rule). Instead, we define the weights using 
ip n in the same way we used w n in fll3[) . More specifically, for every integer j = 1, . . . , n, we define 
the function tpj : [— 1, 1] — >• M via the formula 

ri®= (w) 

with ip n the obvious analogue of w n in <jXTJ) . We observe that, for every integer i,j — 1, . . . , n, 

v 3 iti) = i i (is) 

1 1=3, 
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analogous to (fTU|) . Viewed as a function on the whole real line, each ipj is bandlimited with the 
same band limit c (see, for example, [25], [25], or Theorem 19.3 in [3T]). We define the weights 
W\ , . . . , W n via the formula 

Wj = J_<P&) dt, (16) 

for every j = 1,2, ... ,n (note the analogy with (fT3"|0 . The weights W\, . . . , W n , defined via (fT5)) . 
are different from the solution of the linear system ([?])■ Nevertheless, the resulting quadrature rule 
turns out to satisfy ([5]), provided that e is of order |A n | (see Theorem [TH in Section |4~T1 for a more 
precise statement). 

The analysis of this issue is somewhat long and involved; the reader is referred to [35J, [35] for 
details and proofs. On the other hand, the underlying ideas are relatively simple: Section I2T41 below 
contains a short overview of this analysis. 



2.4 Overview of the Analysis 

The following observation lies at the heart of the analysis: for any band limit c > and any integer 
n > 0, the reciprocal of tp n can be approximated by a rational function with n poles in (—1,1) 
up to an error of order |A„|, where X n is the nth eigenvalue of the integral operator F c (see @ in 
Section[lJ. In other words, the reciprocal of ip n resembles the reciprocal of a polynomial of order n, 
in the following sense. 

If P is a polynomial with n simple roots z\, . . . , z n in (—1, 1), then the function z — > (P(z)) 1 is 
meromorphic in the complex plane; moreover, 

1 ™ 1 

W) ^n*,) ■(*-*>)' (17) 

for all complex z different from z\, . . . ,z n (this is a special case of the well known Cauchy's integral 
formula: see, for example, |31]V Similarly, the function z — > (ip n (z)) is meromorphic; however, it 
has infinitely many poles, all of which are real and simple (see Remark[5]in Section [3. ip . and exactly 
n of which lie in (—1, 1) (see Theorem[T]in Section l3~Tj) . Suppose that the roots of ip n in (—1, 1) are 
denoted by t\ < • • • < t n . It turns out that 

1 ™ 1 

T77T = E TTuYT, H + ° (|A " I) < (18) 

for all real —1 < t < 1 (note the similarity between ([TTf and ([TBI ). In other words, ([T8"]) means 
that the reciprocal of ip n differs from a certain rational function with n poles by a function whose 
magnitude in the interval [—1,1] is of order |A„|. A rigorous version of (fT51) is provided by Theorem^ 
in Section I3TT1 (its proof is somewhat involved; see [2S], for details). More specifically, according 
to this theorem, 



< |A„| ^24 • log ( JL) + 130 • ( X „) 1/4 ) , (19) 



for all real — 1 < t < 1, where Xn is the nth eigenvalue of the prolate differential operator (see 
Theorem [3] in Section I3TTT) . 
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The identity (fTSI is related to the quadrature, discussed in Section 12.31 above, in the following 
way. Multiplying both sides of (fT%| by i() n {t) and using (fl"4"|). we obtain 



1 = ipx{t) + ■■■ + ip n (t) + Mt) ■ O (|A„|) 



(20) 



In other words, ip\, . . . , tp n constitute a partition of unity on the interval [—1, 1], up to an error of 
order |A„|. We integrate both sides of (|2T)|) over [—1,1] and use Theorem [1] in Section I3TT1 and ([H 
in Section |2~31 to obtain 



W 1 + --- + W n = 2 + 0(\\ n \), 

where W\, . . . , W n are the weights of the quadrature rule (see Section |4"31 for more details). 
Suppose now that m ^ n is an integer. We multiply both sides of (|20l) by tp m to obtain 



(21) 



iM*) = ( S^nW-ViW ] + 4>m(t) ■ 4>n(t) ■ O {\\ n \) . 



(22) 



A detailed analysis of a combination of (fl9|) and (1221) leads to the conclusion that, for all integer 

< m < n, 



i) m (t)dt-J2^m(tj)-Wj 



< |An|- ( 24-l0g^ + 6- X n 



(23) 



(see Theorem [TJ] in Section |4~T1 and also [5S], [35] for more details ). 

According to (|2"3"j) . the quadrature rule of order n integrates the first n PSWFs to precision 
of order |A n | (see also © in Section |2~2"|) . It remains to establish for what values of n this error 
is smaller than a predetermined e > 0. Theorem [TH] from Section 14.21 provides an answer to this 
question: namely, if 



n > 



2c 



10 + ^-log(c) + i-logij 'log (|) 



(24) 



then 



3 = 1 



< E, 



(25) 



for all integer < m < n. 

Numerical experiments seem to indicate that the situation is even better in practice: namely, 
to achieve the accuracy e it suffices to pick the minimal n such that |A n | < e, which occurs for 
n « 2c/-7r + 2(logc) • (— loge)/-7r 2 (see Section [7J in particular, Conjectures [3j [4] and Experiment 3 in 
Section I7T|) . 



3 Mathematical and Numerical Preliminaries 



In this section, we introduce notation and summarize several facts to be used in the rest of the 
paper. 
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3.1 Prolate Spheroidal Wave Functions 

In this subsection, we summarize several facts about the PSWFs. Unless stated otherwise, all these 
facts can be found in [37], [3D], [IS], [33J, [T3], [21], US- 
Given a real number c > 0, we define the operator F c : L 2 [— 1, 1] — > L 2 [— 1, 1] via the formula 

F c [ip] (x) = [ v{tY cxt dt. (26) 



Obviously, F c is compact. We denote its eigenvalues by Ao, X\, ■ • ■ , A n , . . . and assume that they 
are ordered such that |A„| > |A n +i| for all natural n > 0. We denote by ip n the eigenfunction 
corresponding to A„. In other words, 



l 

icxi 



X n ip n {x)= J ip n (t)e lcxt dt, (27) 

for all integer n > and all real — 1 < x < 1. We adopt the convention^ that H^rilli 2 ^ 1,1] = 1- The 
following theorem describes the eigenvalues and eigenfunctions of F c . 



Theorem 1. Suppose that c > is a real number, and that the operator F c is defined via (I26p 
above. Then, the eigenfunctions ipo,ipi, ... of F c are purely real, are orthonormal and are complete 
in L 2 [—1, 1]. The even-numbered functions are even, the odd-numbered ones are odd. Each function 
ip n has exactly n simple roots in (—1,1). All eigenvalues X n of F c are non-zero and simple; the 
even-numbered ones are purely real and the odd-numbered ones are purely imaginary; in particular, 
A n = i"|A„|, for every integer n > 0. 

We define the self-adjoint operator Q c : L 2 [— 1, 1] — > L 2 [— 1, 1] via the formula 

^ r w x 1 f 1 sia(c(x — t)) , , , 

Q c [<p] (x) = - ' <P® dt. (28) 

7T l_t X-t 



Clearly, 



Qc M (x) = X[-i,i](x) ■ [x[-cx](0-?M(0] (.x), (29) 



where 3^ : i 2 (R) — > L 2 (M.) is the Fourier transform, and X[-a.a] '■ K — > M is the characteristic function 
of the interval [—a, a], defined via the formula 

/ \ J 1 -a < x < a, 
X[-a,a]W = i„ . (30) 

10 otherwise, 

for all real x. In other words, Q c represents low-passing followed by time-limiting. Q c relates to F c , 
defined via ([2"6"| . by 

Qc = f - • F* ■ F c , (31) 

Z7T 

and the eigenvalues /i„ of Q n satisfy the identity 

= ■ |A„| 2 , (32) 



1 This convention agrees with that of I37| . |30| and differs from that of I33| . 
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for all integer n > 0. Obviously, 

fln<l, (33) 



for all integer n > 0, due to ([29]). Moreover, Q c has the same eigenfunctions ip n as F c . In other 
words, 

, , . If 1 sin (c(x — t)) , , . , . ,. 

^«(a;) = - / ' Mt) dt, (34) 

for all integer n > and all — 1 < x < 1. Also, Q c is closely related to the operator P c : L 2 (U) — > 
i 2 (R), defined via the formula 

„ r i / x 1 /"°° sin (c (x — t)) , . , ,„_. 
P c ^ x = - / vv yy </> (*) dt, 35 

which is a widely known orthogonal projection onto the space of functions of band limit c > on 
the real line K. 

The following theorem can be traced back to [T3] : 

Theorem 2. Suppose that c > and < a < 1 are positive real numbers, and that the operator 
Q c : L 2 [—1, 1] — > L 2 [—1, 1] is defined via (|28p above. Suppose also that the integer N(c,a) is the 
number of the eigenvalues fx n of Q c that are greater than a. In other words, 

N(c, a) = max{fc = 1, 2, . . . : ^ifc_i>Q!}. (36) 

Then, 

N (c, a) = -+(-^\og^-) log c + O (log c). (37) 



According to (l37|) . there are about 2c/ it eigenvalues whose absolute value is close to one, order logc 
eigenvalues that decay rapidly, and the rest of them are very close to zero. 

The eigenfunctions ip n of Q c turn out to be the PSWFs, well known from classical mathematical 
physics |20j . The following theorem, proved in a more general form in |34j . formalizes this statement. 

Theorem 3. For any c > 0, there exists a strictly increasing unbounded sequence of positive numbers 
Xo < Xi < • • • such that, for each integer n > 0, the differential equation 

(1 - x 2 ) ■ ip"(x) ~ 2x ■ %j)'{x) + (xn - c 2 x 2 ) ■ i){x) = (38) 

has a solution that is continuous on [—1,1]. Moreover, all such solutions are constant multiples of 
the eigenfunction %p n of F c , defined via (|26[) above. 

Remark 6. For all real c > and all integer n > 0, (|27[) defines an analytic continuation of ipn 
onto the entire complex plane. All the roots of tp n are simple, real, and symmetric about the origin. 
Moreover, ij) n has infinitely many roots in (l,oo). In addition, the ODE Q38p is satisfied for all 
complex x. 

Many properties of the PSWF ip n depend on whether the eigenvalue \n of the ODE (j38|) is greater 
than or less than c 2 . In the following theorem from [21] , |22j . we describe a simple relationship 
between c, n and Xn- 
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Theorem 4. Suppose that n > 2 is a non-negative integer. 

• If n < (2c/ n) — 1, then Xn < c 2 . 

• If fi > (2c/ it), then Xn > c 2 - 

• If (2c/ it) — 1 < n < (2c/ w), then either inequality is possible. 

In the following theorem, upper and lower bounds on \n in terms of c and n are provided. 
Theorem 5. Suppose that c > is a real number, and n > is an integer. Then, 

n(n + 1) < Xn <n(n + 1) + c 2 . (39) 



It turns out that, for the purposes of this paper, the inequality ([39)) is insufficiently sharp. Tighter 
bounds on Xn are described in the following theorem (see |21j . [22]). 



Theorem 6. Suppose that n > 2 is an integer, and that Xn > c 2 ■ Then, 



2 

n < — 



2+2 



Xn - C-t 
l-t 2 



dt <n + 3. 



(40) 



In the following theorem from |23j , |24j , we provide an upper bound on | A n | in terms of n and c. 
Theorem 7. Suppose that c > is a real number, and that 

c > 22. (41) 
Suppose also that 5 > is a real number, and that 



Suppose, in addition, that n is a positive integer, and that 

2c 2 r , /4e7rc\ 
n>- + - 2 - *-log(— J- 

Suppose furthermore that the real number £(n, c) is defined via the formula 

6 



£ (n, c) — 7056 • c • exp 



-6 1 



2-kc 



(42) 



(43) 



(44) 



The 



|A„| < £(ra,c). 



(45) 



In the following theorem, we provide a recurrence relation between the derivatives of tp n of 
arbitrary order (see Lemma 9.1 in [37]). 
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Theorem 8. Suppose that c > is a real number, and that n > is an integer. Then, 

(1 - t 2 ) - 4t<(t) + ( Xn - c 2 t 2 - 2) yj'Jt) - 2c 2 t^ n (t) = (46) 

for all real t. Moreover, for all integer k > 2 and all real t, 

(1 - t 2 ) V^ +2) (t) - 2 (k + 1) (t) + (X« - fc (fc + 1) - c 2 t 2 ) (t) 

- C 2 kt^~V (t) _ c 2 fc ( fc _ !) (t j = ^ (47) 



The following theorem asserts that, on the interval [—1,1], the difference between the reciprocal 
of if) n and a certain rational function with n poles is of order |A„|. Its proof can be found in [25j . 

M- 

Theorem 9. Suppose that c > 30 is a real number, that n is a positive integer, and that 

2c . , 

n > — + 7. (48) 

Suppose furthermore that —1 < t\ < ■■■ < t n < 1 are the roots of ip n in (—1,1), and that the 
function 5 : [— 1 , 1] — > R is defined via the formula 

1 " 1 

m = AM~^[^)-(t- tj y (49) 

for all real —1 <t< 1. Then, 

W)\ < |An| • (^24 • log ( JL) + 130 • ( X „) 1/4 ) , (50) 

for all real — 1 < t < 1. 

Remark 7. Suppose that the function 5 : [—1, 1] — >• R is defined via (I49|) . 7/n is even, i/ien <5 is an 
even function. If n is odd, then 5 is an odd function. 

3.2 Legendre Polynomials and PSWFs 

In this subsection, we list several well known facts about Legendre polynomials and the relationship 
between Legendre polynomials and PSWFs. All of these facts can be found, for example, in [9], [37] . 

DP- 

The Legendre polynomials Pq, Py, P%, . . . are defined via the formulae 

Po(t) = 1, 

Pi(t) = t, (51) 

and the recurrence relation 

(k + l)P k+1 {t) = (2k + l)tP k (t)-kP k ^(t), (52) 

for all k — 1, 2, . . . . Even Legendre polynomials are even functions, and odd Legendre polynomials 
are odd. The Legendre polynomials {P k }kLo constitute a complete orthogonal system in L? [—1,1]. 
The normalized Legendre polynomials are defined via the formula 



P k (t)=P k (t)-y/k + 1/2, (53) 
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for all k = 0, 1, 2, The L 2 [—1, l]-norm of each normalized Legendre polynomial equals to one, 



i.e. 



f^Pkit)) 2 dt = l. 



(54) 



Therefore, the normalized Legendre polynomials constitute an orthonormal basis for L 2 [— 1 , 1] . In 
particular, for every real c > and every integer n > 0, the prolate spheroidal wave function ?/)„, 
corresponding to the band limit c, can be expanded into the series 



V„(z) = X>*° ' = X>*° ' (55) 

fe=0 fc=0 



for all —1 < x < 1, where /3q™\ /3[ n \ . . . are defined via the formula 



r-l 

(3 { h n) = I il> n {x) -I\{x) dx, (56) 



and . . . are defined via the formula 



i n) = ^ ■ Vfc + 1/2 = (k + 1/2) • / Mx) ■ Pk{x) dx, (57) 



for all k = 0, 1, 2, Due to the combination of Theorem Q] in Section IO with (jS"4"]) . (jS"5"j) . 

2 / , ,\ 2 



(/#°) + (/#°) +"- = l. (58) 



For any integer n > 0, the sequence (3^, /3[ n \ . . . satisfies the recurrence relation 





+ A , 2 


■& ] 


= Xn 






] +A h3 




= Xn 




A ktk . 2 .^ 2 +A k , k -^ + 


Ak,k+2 ■ 




= Xn 


• 4" 



(59) 

for all k = 2, 3, . . . , where A k , k , A k+ 2,k, A k _ k +2 are defined via the formulae 

2fc(fc+l)-l , 
^- fc ( fc + 1 ) + ( 2fc + 3) (2 i-l) - C ' 
A A (fc + 2)(fc + 1) a 

(2/c + 3)V(2fc + l)(2fc + 5) 

for all fc = 0, 1, 2, ... . In other words, the infinite vector (/3q^ , (3^ , . . . J satisfies the identity 

(A- Xn I)-(ti n \i3{ n \...) T = 0, (61) 

where I is the infinite identity matrix, and the non-zero entries of the infinite symmetric matrix A 
are given via (j60l) . 

The matrix A naturally splits into two infinite symmetric tridiagonal matrices, A even and A odd , 
the former consisting of the elements of A with even-indexed rows and columns, and the latter 
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consisting of the elements of A with odd-indexed rows and columns. Moreover, for every pair of 
integers n, k > 0, 



4 n) =0, if k + n is odd, 



(62) 



due to the combination of Theorem Q] in Section I3~T1 and (l56l . In the following theorem (that appears 
in [3 7) in a slightly different form), we summarize certain implications of these observations, that 
lead to numerical algorithms for the evaluation of PSWFs. 

Theorem 10. Suppose that c > is a real number, and that the infinite tridiagonal symmetric 
matrices A even and A odd are defined, respectively, via 



and 



A 



odd 



( A),0 Aq_2 
Alfi A%,2 ^2,4 

^4,2 ^4,4 ^4,6 



V 



(A h i Ai,3 
A^,l Ay,^ ^3 : 5 

^5,3 ^5,5 A^j 



■J 



(63) 



V 



(64) 
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where the entries Akj are defined via (|60|) . Suppose also that the infinite vectors Peven € I 2 and 
^odd ^ ^ 2 are defined, respectively, via the formulae 

p£L = ■ M n) , • ■ • ) T , tdl = (M n) ■ ^ (65) 

where /3q , /3{ , . . . are defined via (|56[) . If n is even, then 

A ' ft even ~ Xn ' ft even- (66) 

If n is odd, then 



A odd , o(n) _ o(n) 
71 Podd — Xn P odd - 



(67) 



Remark 8. We define the infinite vector (3^ € Z 2 £o be equal to Peven> if n * s even, or to P^ d , 
if n is odd. In this notation, fj(°\ fj( 2 \ . . . are the eigenvectors of A even , and fj^,/3^\... are the 
eigenvectors of A odd . 

Remark 9. While the matrices A even and A odd are infinite, and their entries do not decay with 
increasing row or column number, the coordinates of each eigenvector (3^ decay superexponentially 
fast (see e.g. Jff7| / for estimates of this decay). In particular, suppose that we need to evaluate 
the first n + 1 eigenvalues Xo> • • • > Xn and the corresponding eigenvectors ■ ■ ■ , j3( n ' numerically. 
Then, we can replace the matrices A even , A odd in (|66[) . (1671) . respectively, with their N x N upper left 
square submatrices, where N is of order max{n,c}, and solve the resulting symmetric tridiagonal 
eigenproblem by any standard technique (see, for example, \36f . see also \31^ for more details 
about this numerical algorithm). The CPU cost of this procedure is 0(n 2 ) operations. 
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The Legendre functions of the second kind Qo, Qi, Q2, 



are defined via the formulae 



Qo(t) = - lo S 

Q 1 (t) = ^log^-l, (68) 

and the recurrence relation 

(k + 1) Qfc+i(t) = (2k + 1) tQ k (t) - fcQ fc _i(t), (69) 

for all fc = 1,2, ... . We observe that the recurrence relation is the same as the recurrence 
relation (|52j). satisfied by the Legendre polynomials. In addition, for every integer k — 0,1,2,..., 
the fcth Legendre polynomial and the feth Legendre function of the second kind Qk are two 
independent solutions of the second order differential equation 

(1 - t 2 ) ■ y"(t) - 2t ■ y'(t) + k(k + 1) ■ y(t) = 0. (70) 



Remark 10. Suppose that —1 < x < 1 is a real number, and that n > is an integer. Com- 
bining (f5"Tj) . (|52p . (JHSJ) 7 (1551 gives a numerical procedure for the evaluation of Pq(x), ... ,P n (x) and 
Qq(x), . . . ,Q n (x) to high precision. This procedure is stable, and requires 0{n) operations (see, for 
example, for more details). 



3.3 Priifer Transformations 

The classical Priifer transformation of a second-order ODE is a well known analytical tool for the 
study of the oscillatory properties of its solutions (see, for example, [IH], [()])■ Recently, a minor 
modification of Priifer transformation was demonstrated to be also a convenient numerical tool (see 
[8]). In the following theorem, we summarize several properties of this transformation, applied to 
the prolate ODE (38J) (see 0, HB, HU for details). 

Theorem 11. Suppose that n > 2 is an integer, and that \n > c 2 . Suppose also that the functions 
f,v : (—1, 1) — > M are defined, respectively, via the formulae 

and 

1 ( t c 2 t \ 

for all real —1 < t < 1. Suppose furthermore that t\ is the minimal root of tpn in (—1, 1), and that 
the function : (— 1, 1) — > R is the solution of the differential equation 

6'{t) = f(t) + v(t) ■ sin(20(i)) (73) 



with the initial condition 

0(h) - \ (74) 

Then, has the following properties: 
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extends continuously to the interval [—1,1], and, moreover, 

0(-l)=O, (75) 



7171 



0(0) = -y, (76) 
6»(1) = 7m. (77) 

For any real — 1 < t < 1 suc/i i/iai ipn(t) =/= 0, 



6(t) = atan(-J 1 ■ ^§ ) + m(t) ■ tt, (78) 



X„ - c 2 t 2 ip n (t) 

where m(t) is the number of the roots of if> n in the interval (— l,t). 

• For each integer i = 1, . . . , n, 

6(ti) = (i-±y*, (79) 

where ti,...,t n are the roots of ip n in (—1, 1). 

• For all real — 1 < t < 1, 

6'{t) > 0. (80) 
In other words, is monotonically increasing. 

The following theorem is closely related to Theorem [TT] (see [21], [22] for more details). 

Theorem 12. Suppose that the function 8 : [£i,£„] — > K that of Theorem ] 111 Suppose also that the 
function s : [7r/2,7r • (n — 1/2)] — > [t\,t n ] is the inverse of 9. Then, s is well defined, monotonically 
increasing and continuously differentiable. Moreover, for all real tt/2 < r/ < it ■ (n — 1/2), 

s ' {v) = fiM^ + vism-Miv)' (81) 

where the functions f, v are defined, respectively, via (|71[) , (|72p . In addition, for every integer 
i = l,...,n, 

« |) " 7r ) =**' (82) 

and also 

.(=)-a (83, 



3.4 Numerical Tools 

In this subsection, we summarize several numerical techniques to be used in this paper. 
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3.4.1 Newton's Method 

Newton's method solves the equation f(x) = iteratively given an initial approximation xq to the 
root x. The nth iteration is defined by 

Xn X n —i — — - . (^4) 

The convergence is quadratic provided that J is a simple root and xq is sufficiently close to i. More 
details can be found e.g. in [5]. 

3.4.2 The Taylor Series Method for the Solution of ODEs 

The Taylor series method for the solution of a linear second order differential equation is based on 
the Taylor formula 

is) { \ 

u{x + h) = Y, ! —f l W +0(h k+1 ). (85) 

3=0 J ' 

This method evaluates u{x + h) and u'(x + h) by using (|8"5|) and depends on the ability to compute 
u^\x) for j = 0, ... , k. When the latter satisfy a simple recurrence relation such as (|4"T)) and hence 
can be computed in 0(k) operations, this method is particularly useful. The reader is referred to 
[5] for further details. 

3.4.3 A Second Order Runge-Kutta Method 

A standard second order Runge-Kutta Method (see, for example, [5]) solves the initial value problem 

y(to) = Vo, y'(t) = f(t,y) (86) 

on the interval to < t < io + L via the formulae 

ti+i = £» + h, 

h+i = hf (t i+1 ,yi + ki) , 

y i+ i=yi + {ki + k i+ i)/2 (87) 
with i = 0, . . . , n, where h and fco are defined via the formulae 

h=- k = f(t ,y ). (88) 
n 

This procedure requires exactly n + 1 evaluations of /. The global truncation error is 0(h 2 ). 

3.4.4 Shifted Inverse Power Method 

Suppose that n > is an integer, and that A is an n by n real symmetric matrix. Suppose also that 
01 < 02 < • • • < o>i are the eigenvalues of A. The Shifted Inverse Power Method iteratively finds 
the eigenvalue Cfe and the corresponding eigenvector Vk G R™, provided that an approximation A to 
ak is given, and that 

I A — o-fc I <max{|A-<7 i | : j ^ k} . (89) 
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Each Shifted Inverse Power iteration solves the linear system 

(A -\ j I)-x = Wj (90) 

in the unknown x € R™, where Xj and Wj € R™ are the approximations to a k and v k , respectively, 
after j iterations; the number Xj is usually referred to as "shift". The approximations Xj+i and 
Wj+i G R™ (to (Tfc and Vk, respectively) are evaluated from x via the formulae 

Wj+l = n-Ty = wj +1 ■ A ■ W j+1 (91) 

(see, for example, [5], [35] for more details). 

Remark 11. For symmetric matrices, the Shifted Inverse Power Method converges cubically in the 
vicinity of the solution. In particular, if the matrix A is tridiagonal, and the initial approximation 
X is sufficiently close to a k , the Shifted Inverse Power Method evaluates o~k and Vk essentially to 
machine precision e in O {— log log e) iterations, and each iteration requires 0(n) operations (see e.g 



3.4.5 Sturm Bisection 

In this subsection, we describe a well known algorithm for the evaluation of a single eigenvalue of 
a real symmetric tridiagonal matrix. This algorithm is based on the following theorem that can be 
found, for example, in |36) . [2J. 

Theorem 13 (Sturm sequence). Suppose that n > is an integer, that 



C 



(a 1 
b 2 







b 2 
a-2 




h 



bn-l 









a n -i b n 
b n a n J 



(92) 



is an n by n symmetric tridiagonal matrix, and that none of numbers b2,---,b n is equal to zero. 
Suppose also that the polynomials p^i,po, . . . ,p n are defined via the formulae 

P-i(x) = 0, p Q (x) = 1 

and 

p k (x) = (a k - x)p k -i{x) - b 2 k p k - 2 (x), 



(93) 



(94) 



for all real x and every integer k = 2, . . . , n. Suppose furthermore that a is a real number, and that 
the integer A{o~) is defined as the number of positive elements in the finite sequence 



Poio-)pi(o-), pi(a)p 2 (o-), p n -i(o-)p n (a). 
Then, the number of eigenvalues of C that are strictly larger than a is precisely A(a). 



(95) 



Remark 12. Suppose now that n > is an integer, and C is an n x n real symmetric tridiagonal 
matrix, such as (|92[) . Theorem ] ISA yields a numerical scheme for the evaluation of the kth smallest 
eigenvalue a k of C . This scheme is known in the literature as "Sturm Bisection" . Provided that two 
real numbers xq and yg are given such that 



X < CTfc < yo, 



(96) 
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Sturm Bisection requires 



operations to evaluate o~k to machine precision (see, for example, /<%] /. fBj for more details). 

4 Analytical Apparatus 

The purpose of this section is to provide the analytical apparatus to be used in the rest of the paper. 
More specifically, we define a PSWF-based quadrature rule and list several of its properties. 

The principal result of this section is Theorem \W\ The reader is referred to [35] , [33] for the 
detailed analysis of all the tools listed in this section. 

Throughout this section, the band limit c > is assumed to be a positive real number. Also, 
for any integer n > 0, we denote by ip n the nth PSWF corresponding to the band limit c (see 
Section O). 



Definition 2. Suppose that n > is an integer, and that 

-1< h < t 2 < ■ ■ ■ < t n < 1 (98) 

are the roots of ip n in the interval (—1,1). For each integer j = l,...,n, we define the function 
ipj : [— 1, 1] — > M via the formula 

In addition, for each integer j = 1, . . . ,n, we define the real number Wj via the formula 

w >--l>^--^kL^ 

We refer to the pair of finite sequences 

S n = (ti,...,t n ,Wi,...,W n ) (101) 

as the "PSWF-based quadrature rule of order n". The points t\, . . . , t n are referred to as the quadra- 
ture nodes, and the numbers W\, . . . ,W n are referred to as the quadrature weights (see (|3" j) . (pT | in 
Section^). We use S n to approximate the integral of a bandlimited function f over the interval 
[—1,1] by a finite sum; more specifically, 

f f{t)dt^Y, w r f(h)- (102) 

We refer to the number S n (f) defined via the formula 



8n(f) 

as the "quadrature error". 



/I n 
_1 i=i 



(103) 
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4.1 Quadrature Error and its Relation to |A n | 

Suppose now that n is a positive integer, and that / : [—1,1] — > C is an arbitrary bandlimited 
function (with band limit c). Suppose also that S n is the PSWF-based quadrature rule of order n 
(see l|101|) in Definition [2]) . One of the principal goals of this paper is to investigate the quadrature 
error 5 n (f) defined via (|103[) . The reader is referred to Section [7] for the results of several related 
numerical experiments. 

The following theorem, illustrated in Table [TJ provides an upper bound on 8 n {ip m ), for any 
integer m = 0, . . . , n — 1. This theorem is illustrated in Tableland in Figure [5] (see Experiment 2 
in Section fTTTj) : see also Conjecture [3] and Remark l37l in Section I77L1 



Theorem 14. Suppose that c is a positive real number, and that 

c > 30. (104) 

Suppose also that n > and < m < n — 1 are integers, and that 

2c 

n> — + 5. (105) 
Suppose further that 5 n (?p m ) is defined via (|103p . Then, 



/l n 
1p m (s) ds-^ Wj ■ 1p m (tj) 
- 1 .7=1 



: |A„|- (24.|,„( J-,) (i . x „). (],„;, 



where \ n ,Xn are those of (|27p . (|38p in Section \S.l\ respectively. 

4.2 Quadrature Error and its Relation to n and c 

In Theorem [TH we established an upper bound on the quadrature error 8 n (i/j m ) (see (I103p and (|106p 
in Theorem 1 14p . However, this bound depends on Xn an d A„. In particular, it is not obvious how 
large n should be to make sure that the quadrature error does not exceed a prescribed e > 0. In 
this subsection, we eliminate this inconvenience. 

The following theorem is illustrated in Table [H (see Experiment 3 in Section PTTj) . 

Theorem 15. Suppose that c,e are positive real numbers such that 

c> 30 (107) 

and 

< log - < -^J= ■ c - 3 • log(c) - log(6 5 ■ 14340). (108) 
Suppose also that the real numbers a, v(a) are defined via the formulae 

a = • flog ~ + 3 • log(c) + log(6 5 • 14340) j (109) 



, . 2c a , ( 16ec\ 
"(«) = - + 77Z ' — . ( 110 ) 



f 2i V OL 
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respectively. Suppose furthermore that n > and < m < n — 1 are integers such that 

n>v(a), (111) 

and that 8 n {ip m ) is defined via (|103|) . Then, 



8 n {ip-n 



/l n 
tp m {s) ds-^ ^mit^Wj 
- 1 7 = 1 



< e. 



(112) 



The following theorem is a direct consequence of Theorem [T5] This theorem is one of the 
principal results of the paper. It is illustrated in Table 2] (see Experiment 3 in Section [77T|) . See also 
Conjecture [3] in Section 17.11 



Theorem 16. Suppose that c,e are positive real numbers such that 

c> 60 

and 

< e < 1. 

Suppose also that n > and < m < n are integers, and that 

2c / 3 , , x 1 , 1\ , /c\ 
fl >_ + ^10 + -.Iog( C ) + -.log-J.log(-). 

Suppose furthermore that 8 n {tp m ) is defined via (|103l) in Definition® Then, 



Sn{tprn) = 



/I n 
1p m (s) ds-^ i)m(tj)Wj 
- 1 .7 = 1 



< e. 



(113) 



(114) 



(115) 



(116) 



4.3 Quadrature Weights 

In this subsection, we analyze the weights of the quadrature rule S n (see (|100l) , (|101j) in Section 0} . 
This analysis has two principal purposes. On the one hand, it provides the basis for a fast algorithm 
for the evaluation of the weights. On the other hand, it provides an explanation of some empirically 
observed properties of the weights. 

The results of this subsection are illustrated in Table [5] and in Figure [6] (see Experiment 4 in 
Section £2Jl. 

The following theorem is instrumental for the evaluation of the quadrature weights W± , . . . , W n 
(see ifTUD)) in Definition 0). 

Theorem 17. Suppose that n > is an integer, and that the function <!>„ : (—1,1) — > R is defined 
via the formula 

oo 

*„(*) = 5>jp-Q*(t), (117) 

fc=0 
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where Qk(t) and ctu" are defined, respectively, via (|68[) . (|69p and (|57p in Section \ 3.S\ (compare to 
(|55p m Section \3.!ty) . Then, for every integer j = 1, . . . , n, 

^ = -^py^ = ->-Wr <118) 

where t\, . . . , t n and W\, . . . , W n are, respectively, the nodes and weights of the quadrature rule S n 
in Definition^ 

Theorem[T7]is illustrated in Table[5J We observe that TheoremQTT] describes a connection between 
the weights W\, . . . , W n and the values of $ n at t\, . . . , t n , where the function is defined via (|117p . 

The following theorem states that $„ satisfies a certain second-order non-homogeneous ODE, 
closely related to the prolate ODE ((38)) in Section 137X1 In particular, a recurrence relation between 
the derivatives of <t n of arbitrary order is established (compare to Theorem [5] in Section 13.11) . 

Theorem 18. Suppose that n > is an integer, and that the function $„ : (—1, 1) — >■ M is defined 
via (I117p . Suppose also that the real numbers a^\a^ are defined via (|57p in Section [3~M Then, 

(1 - t 2 ) ■ - 2t ■ K(t) + (Xn ~ c 2 t 2 ) • i n (t) = -c 2 (4 n) t + 4 l) /3) , (119) 

/or aH reaZ —1 < t < 1. ^4feo, 

(1 - t 2 ) ■ $'i'(t) - At ■ K(t) + ( X n - c 2 t 2 - 2) • - 2c 2 i • * n (t) - -c 2 o>\ (120) 

/or a/Z reaZ — 1 < t < 1. Finally, 

(1-t 2 ) $(f+ 2 ) (t) - 2 (fc + 1) i*^ 1 ) (i) + ( Xn - fc (fc + 1) - c 2 t 2 ) * W (t) 

- c 2 kt¥*-v (t) - c 2 k (k - i) ¥ n k - 2) (t) = 0, (121) 

for every integer k > 2 and a?/ reaZ — 1 < t < 1 (compare to (1471) in Section \'3.1\) . 

In the following theorem, we establish the positivity of the weights of the quadrature rule S n in 
Definition [5J 

Theorem 19. Suppose that c is a positive real number, and that 

c > 30. (122) 
Suppose also that n is a positive integer, and that 

n > ^ + 5 ■ log(c) • log (|) . (123) 

Suppose further that W\, . . . , W n are defined via (|100p . Then, for all integer j = 1, . . . , n, 

Wj > 0. (124) 



Remark 13. Extensive numerical experiments (see e.g. Tableland Figure^) seem to indicate that 
the assumption (|123p is unnecessary. In other words, the weights W± , . . . , W n are always positive, 
even for small values of n (at the present time we do not have the proof of this fact). 
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Remark 14. It was observed in J25f . \26$ that, if 1 < j,k < n are integers, then 

(Cfe)) 2 • (1 - ^) ■ Wi = «(*fc)) 2 • (1 - tl) -W k +0 (|A„|) (125) 



(see also Experiment 4 in Section \7.2\ ). We observe that as c — > the quadrature rule in Definition^ 
converges to the well known Gaussian quadrature rule, whose nodes are the roots t\ , . . . , t n of the 
Legendre polynomial P n ( see Section \3.2)) . and whose weights are defined via the formula 

Wj = r, (126) 

] (P^fe)) 2 .(l-^)' 

for every j = 1, .. . ,n (see e.g. Ji/ ; Section 25.4). Thus, ()125[) is not surprising. 

5 Numerical Algorithms 

In this section, we describe several numerical algorithms for the evaluation of the PSWFs, certain 
related quantities, and the quadrature rules defined in Section 0] Throughout this section, the band 
limit c > is a real number, and the prolate index n > is a non-negative integer. 



5.1 Evaluation of Xn and if) n (x), i)' n {x) for — 1 < x < 1 

The use of the expansion of ij) n into a Legendre series (see ((55)) in Section I3.2|) for the evaluation 
of ip n in the interval [—1,1] goes back at least to the classical Bouwkamp algorithm (see More 
specifically, the coefficients /3q™^ , f3[ 7 ^ , . . . of the Legendre expansion are precomputed first (see (|5"oT) , 
(f5"7]) in Section 1572")) . These coefficients decay superalgebraically; in particular, relatively few terms 
of the infinite sum (|55p are required to evaluate ip n to essentially machine precision (see Section [3.21 
in particular Theorem [TU] and Remark |H1 and also [37] for more details). 



5.1.1 Evaluation of Xn and /3 { ' , /3\ . . . 

Suppose now that n > 0, and one is interested in evaluating the coefficients /3q , /3[ m ^ , ... in (|55p . 
for every integer < m < n. This can be achieved by solving two N x N symmetric tridiagonal 
eigenproblems, where N is of order n (see Theorem [TU] and Remark in Section l3~2l and also [37] 
for more details about this algorithm). In addition, this algorithm evaluates xot ■ ■ iXn- Once this 
precomputation is done, for every integer < m < n and for every real — 1 < x < 1 one can evaluate 
ip m (x) m 0(n) operations, by computing the sum (|55D (see, however, Remark |2"T1 below) . 

Suppose, on the other hand, that we are interested in a single PSWF only (as opposed to all the 
first n PSWFs). Obviously, we can use the algorithm above; however, its cost is 0(n 2 ) operations 
(see Remark [pj in Section 13. 2p . In the rest of this subsection, we describe a procedure for the 
evaluation of /?o™\ P\ , ■ ■ ■ and Xn, whose cost is 0(n + clog(c)) operations. 

This algorithm is also based on ThcorcmfTUlin Section l372"1 It consists of two principal steps. First, 
we compute a low-accuracy approximation Xn of Xm by means of Sturm Bisection (see Section l3.4.5l 
(lUBl . (JSTJ and Remark[UJin Section [3721 and also [5]). Second, we compute Xn an d (see ([rJ5]) and 
Remark [3] in Section 1372]) by means of the Shifted Inverse Power Method (see Section l3.4.4l and also 
[36] . [5]). The Shifted Inverse Power Method requires an initial approximation to the eigenvalue; for 
this purpose we use Xn- 

Below is a more detailed description of these two steps. 
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Step 1 (initial approximation \ n of Xn)- Suppose that the infinite symmetric tridiagonal 
matrices A even and A odd are defined, respectively, via (|63p . (|64l) in Section [XU Suppose also that 

is the N x N upper left square submatrix of A even , if n is even, or of A odd , if n is odd. 
Comment. N is an integer of order n (see Remark [HI in Section 1X2"]) . The choice 

AT= 1.1- c + n + 1000 (127) 

is sufficient for all practical purposes. 

• use Theorems 01 [BJ and [BJ in Section 13.11 to choose real numbers xo < yo such that 

XQ < Xn < yo- (128) 

Comment. For a more detailed discussion of lower and upper bounds on Xn, see, for example, 
[2"T) . [2"2] . See also Remark [TBI below. 

• use Sturm Bisection (see Section [3.4.5p with initial values xq, Ho to compute % n . On each step 
of Sturm Bisection, the Sturm sequence (see (IMj) in Theorem [T5)l is computed based on the 
matrix A 1 - 71 * 1 (see above). 

Comment. In principle, Sturm Bisection can be used to evaluate Xn to machine precision. 
However, the convergence rate of Sturm Bisection is linear, and each iteration requires order 
n operations (see Remark [12] in Section I3.4.5[) . On the other hand, the convergence rate of 
the Shifted Inverse Power Method is cubic in the vicinity of the solution, while each iteration 
requires also order n operations (see Remark ITT1 in Section l3.4.4[) . Thus, we use Sturm Bisection 
to compute a low-order approximation x~n to Xn , and then refine it by the Shifted Inverse Power 
Method to obtain x-n to machine precision. 

Remark 15. The use of Sturm Bisection as a tool to compute the eigenvalues of a symmetric 
tridiagonal matrix goes back at least to Jl|/; in the context of PSWFs, it appears in J10\j . 

The cost analysis of Step 1 relies on the following observation based on Theorems [31 [H [SJ [BJ in 
Section 13.11 

Observation 1. Suppose that n > is an integer. 
If < n < 2c/tt, then 

Xn+l~Xn = 0(c). (129) 

If n > 2c/ it, then 

Xn+l ~ Xn = 0(n). (130) 

Remark 16. Due to Theorems^ [5| in Section \3.1\ the inequality 

n-(n + l)< X n< c 2 (131) 

holds for any real c > and all integer < n < 2c/n. In this case, Step 1 requires 0(c ■ log(c)) 
operations, due to the combination of (|129p , (|131[) and Remark in Section \3.4-5\ On the other 
hand, if n > 2c/ V, then the cost of Step 1 is 0(n) operations, due to the combination of Theorems^ 
Remark\Min Section \3~^~5\ and (lT5Bj) . 
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Step 2 (evaluation of \ n and /3^ n '). Suppose now that \ n is an approximation to Xn evaluated in 
Step 1. Suppose also that the integer N is defined via (|127l) above (see also Remark[5]in Section l3~2"j) . 

• generate a pseudorandom vector $ £ M. N of unit length. 

Comment. We use x„ and (3 as initial approximations to the eigenvalue \n an d the corre- 
sponding eigenvector, respectively, for the Shifted Inverse Power Method (see Section f3. 4. 41) . 

• conduct Shifted Inverse Power Method iterations until \n is evaluated to machine precision. 
The corresponding eigenvector of unit length is denoted by j3^ n \ 

Comment. Each Shifted Inverse Power iteration costs O(N) operations, and essentially 0(1) 
iterations are required (see Remark ITT1 in Section l3.4.4l for more details). In practice, in double 
precision calculations the number of iterations is usually between three and five. 

Remark 17. Clearly, the cost of Step 2 is 0(n) operations (see Remark in Section \3.2\ and 
Remark ] 11\ in Section \3.4-4^ - 

Remark 18. Suppose that the coordinates of the vector ^ n > £ R N are defined via (|65[) (see also 
Remark^ in Section \3.2\) . Then, (3^ (evaluated in Step 2 above) approximates /3^ n ' to essentially 
machine precision ( this is a well known property of the Inverse Power Method; see Section \3.4-4\ 
and also J36jj . J5j/ for more details). In other words, 

11/3^-/3^11 <£• ||/3 (ll) || =£, (132) 

where e is the machine accuracy (e.g. e ~ 1D-16 for double precision calculations). In addition, the 
eigenvalue Xn is also evaluated to relative accuracy e. 

5.1.2 Evaluation of ip n (x), ip' n ( x ) f° r — 1 — x — 1? given Xn and /3q™\ /3^ , . . . 

Suppose now that Xn and the coefficients (3^ , /3^ , . . . defined via (|5"o]) have already been evaluated. 
Suppose also that the integer N is defined via (|127p above. 
For any real — 1 < x < 1, we evaluate ipn( x ) via the formula 

2JV 2N 

* n {x) = E P ^ x ) ■ a k° = E P ^ x ) ■ Pk ] ■ Vk + T/2 (133) 

k=0 k=Q 

(compare to (1551 in Section l3~2"j) . Also, we evaluate ip' n (x) via the formula 

2N 2N 

= E PtW ■ 4*° = E Kb) ■ Pk ] ■ VkTl/2. (134) 

k=l k=Q 

Remark 19. Due to the combination of Remark\§\in Section \3.2\ and Remark ] 18\ above, both ip n (x) 
and ip' n ( x ) are evaluated via (I133p . (|134j) essentially to machine precision, for any real — 1 < x < 1 
(also see J37V for more details). 

Remark 20. Due to Remarks ] 161 \17\ above, the cost of the evaluation of Xn and fto^t Pi > . . . via 
Steps 1,2 is 0(n + clogc) operations. Once this precomputation has been performed, the cost of 
each subsequent evaluation ofip n (x), ip' n {x), for any real —1 < x < 1, is 0(n) operations, according 
to (fT33)l . (fT34|) and RemarkXMin SectionlSTE 

Remark 21. Once Xn and (3^, f3[ n \ . . . have been evaluated, one does not have to use (|133[) , 
(|134l) . to compute tf) n (x), ip' n ( x ) an arbitrary point x in [—1,1]. Instead, the cost of evaluating, 
say, ip n (x) can be brought down from 0{n) to 0(1) (see Remark\£fJ\ in Section [573]) . 
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5.2 Evaluation of A n 

Suppose now that n > is an integer, and that one needs to evaluate the eigenvalue A„ of the 
integral operator F c (see ([2"6"]) in Section 13.11) . Due to the combination of ([2"r?]) and Theorem [T] in 
Section l3~Tl if n is even, then ip n (0) ^ 0, and 



1 



K = —j^ / Mt) dt; (135) 



for odd n, 



C(0) 



A„ = — — / t ■ ip n (t) dt. (136) 



i 



The formulae (|135l) and (|136[) provide an obvious way to calculate A„ for even and odd n, respectively, 
via numerical integration. In fact, when |A„| is relatively large, such procedure is quite satisfactory. 
More specifically, if n < 2c/ V, then |A„| « y/2Tr/c, and A„ can be calculated via (|135[) . (| 136[) to high 
relative precision (see Theorems [3] [7J in Section 13.11 and Remark [T5] in Section 15.11 see also [37] for 
more details). On the other hand, we observe that HV'nllL 2 !-!,!] = 1, due to Theorem[T]in Section I3TT1 
As a result, when |A„| is small, the formulae (I135[) . (|136p are unsuitable for the evaluation of A„ 
via numerical integration, due to catastrophic cancellation. For example, if |A„| < e, where e is the 
machine precision, the formulae f|135[) . (|136[) produce no correct digits at all. 

The standard way to overcome this obstacle for numerical evaluation of small A^s is to calculate 
all the ratios Ao/Ai, . . . , A n /A n _i (see, for example, [H], [33], [31]); this turns out to be a well- 
conditioned numerical procedure (see [37] for more details). Then, Ao is evaluated via (|135l) above, 
and the eigenvalues Ai, . . . , A„ are evaluated via the formula 

A m =A -^ (137) 

for every integer m = 1, . . . , n. 

Suppose, on the other hand, that one is interested in a single A„ only (as opposed to all the 
first n eigenvalues). Obviously, A„ can be evaluated via (I137P from the ratios as described 

above; however, it requires at least 0(n 2 ) operations (see [37]). 

Unexpectedly, it turns out that X n can be obtained to high relative accuracy in O(l) operations 
as a by-product of the algorithm described in Section 15.11 More specifically, suppose that the 
coefficients ^ n \/3[ n) ,. . . are defined via We combine (|T35l) . (fT36) above with ([2T ]) . ([ST]) . (j53 |) . 
([5T>]) . ([57]) to make the following observation. 

Observation 1. If n is even, then ip n (0) ^ 0, and 

A " = w/>" (t)<i< = W' <138) 



1 ^2 icp[ n) 



If n is odd, then ip' n (0) ^ 0, and 

K =mU' M) * = ^rm- (139) 

Remark 22. Obviously, the cost of evaluating A„ from ip n (0), /3g™^ via (|138|) (7or even or from 
ip' n (0), /3± via (|139p (/or odd n) is 0(1) operations. 

Remark 23. £>ue Remarks ] 2(A and (|138p . p39p . a single X n can be evaluated as a by-product 
of the procedure described in Section \5.1[ at the total cost of O (n + clog(c)) operations. 
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Remarks [221 [231 describe the cost of the evaluation of A„ via (|138[) , (|139[) . To describe the accuracy 
of this procedure, we start with the following observation. 

Observation 2. Due to Remark \T^\ X n is evaluated to the same relative accuracy as (3^ (for 
even n) or as f3± (for odd n). According to Q132p in Remark [TS1 the algorithm of Section [5TT1 
evaluates the vector /3' n ) to relative accuracy e, where e is the machine precision. However, this 
means that a single coordinate of f}( n > is only guaranteed to be evaluated to absolute accuracy e. 
More specifically, the inequality 



(■») 



i n) 



< 



(140) 



holds for every integer k = 0, . . . , N, where N is defined via (|127[) in Section 15.11 and $^ is the 
numerical approximation to /3i . In general, the inequality (|140[) can be rather tight; as a result, 



if, for example, |/3q | < £ /10j then, apriori, we cannot expect /3q" - ' to approximate (3q L> to any digit 
at all! 

In practical computations, it is sometimes desirable to evaluate extremely small A„'s (e.g. |A n | f=s 
1D-50). Observation 2 seems to suggest that, in such cases, the evaluation of A„ via the procedure 
described above is futile due to disastrous loss of accuracy. 

Fortunately, it turns out that the algorithm described in Section [5TT1 always evaluates (3^ , (3^ to 
high relative accuracy, regardless of how small they are. This is a consequence of a more general (and 
somewhat surprising!) phenomenon studied in detail in [27], |28j . We summarize the corresponding 
results in the following theorem. 

Theorem 20. For a certain class of real symmetric tridiagonal matrices, the coordinates of their 
eigenvectors are defined to high relative precision. Moreover, the matrices A even , A odd defined, re- 
spectively, via (|63p . (|64l) in Section ] 3. 2[ belong to this class. 

In the following theorem, we summarize implications of Theorem[2n]for the evaluation of (3^ , 
via the algorithm in Section 15.11 (the proof of a slightly modified version of this theorem appears in 

[23, [25]). 

(n) (n) 

Theorem 21. Suppose that c > is a real number, that n > is an integer, and that pg , p\ are 
defined via ((55)) in Section [?TM Then, the algorithm described in Section [5A\ evaluates (3^,(3^ to 
high relative accuracy. More specifically, 



?(") 



?(n) 



/t 3 - 



in) 



for even n, and 



,(») 



0\ 



(n) 



< 10 •£ 



< 10-e-c 



(141) 



(142) 



for odd n, where (3^,(3^ 
machine accuracy (e.g. e f= 



?(n) 



respectively, and e is the 



are the numerical approximation to 0q 
1D-16 for double precision calculations) . 

Remark 24. The algorithm described in Section \5.1\ evaluates the eigenvectors f3^ n ' by the Shifted 
Inverse Power Method (see Section \3.4-4^ - H turns out that the choice of method is important in this 
situation: if, for example, these eigenvectors are evaluated via the standard and well known Jacobi 
Rotations (rather than Inverse Power), the small coordinates exhibit the loss of accuracy expected 
from (|140l) (see \2T$ , \28\/ for more details about this and related issues). 
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Remark 25. Due to the combination of Remark ] 19\ in Section \5. 1[ Observation 2 above and The- 
orem \21l the algorithm of this section evaluates X n to high relative accuracy. More specifically, at 
most 1 + log 10 (c) decimal digits are lost in the evaluation of X n . 

5.3 Evaluation of the Quadrature Nodes 

Suppose that n > is an integer, and that the quadrature rule S n is defined via fllOlf) in Section |U 
According to the nodes of S n are precisely the n roots t\, ■ . ■ , t n of ip n in the interval (— 1, 1). 

In this section, we describe a numerical procedure for the evaluation of these quadrature nodes. 
This procedure is based on the fast algorithm for the calculation of the roots of special functions 
described in 0. It combines Prufer's transformation (see Section l33|) . Runge-Kutta method (see Sec- 
tion l3.4.3"]) and Taylor's method (see Section l3".4.2l) . This algorithm also evaluates tp' n (t\), . . . , tjj' n (t n ). 
It requires 0(n) operations to compute all roots of ip n in (— 1, 1) as well as the derivative of ip n at 
these roots. 

A short outline of the principal steps of the algorithm is provided below. For a more detailed 
description of the algorithm and its properties, the reader is referred to [5]. 
Suppose that i m ; n is the minimal root of ip n in [0, 1). 

Step 1 (evaluation of t m i n ). If n is odd, then 

train = t(n+l)/2 = 0, (143) 

due to Theorem Q] in Section 13.11 On the other hand, if n is even, then 

train = *(n+2)/a > 0- (144) 

To compute t m [ n in the case of even n, we numerically solve the ODE (|5T|) with the initial condition 
(1531 in the interval [7rn/2,7r- (n + l)/2], by using 20 steps of Runge-Kutta method described in 
Section 13.4.31 The rightmost value t m i n of the solution is a low-order approximation of t m i n = 
t(n+2)/2 (see (|82|). (j!44jl ). Then, we evaluate t m i n to machine precision via Newton's method (see 
Section [3.4.ip . using i m i n as an initial approximation to t m i n . On each Newton iteration, we evaluate 
ip n and ip' n by using the algorithm of Section lOl fsee ()133[) . (|134l) ). 

Observation 1. The point i m i n approximates t m i n to at least three decimal digits (see Sec- 
tion 13. 4. 3D . Since Newton's method converges quadratically in the vicinity of the solution, only 
several Newton iterations are required to obtain i m j n from i m j n to essentially machine precision (see 
[8] for more details). In our experience, the number of Newton iterations in this step never exceeds 
four in double precision calculations (and never exceeds six in extended precision calculations). We 
combine this observation with Remark [201 in Section 15.11 to conclude that the total cost of Step 1 is 
0(n) operations. 

Step 2 (evaluation of ip' n (t m i n )). We evaluate ^(t m in) to machine precision via (|134[) in Sec- 
tion o 

Observation 2. Due to RemarkJ^Hlin Section lixll the cost of Step 2 is 0(n) operations. 
The remaining roots of ip n in (i m i n , 1) are computed one by one, as follows. Suppose that 
n/2 < j < n is an integer, and both tj and ip' n {tj) have already been evaluated. 

Step 3 (evaluation of tj + i and ip' n (tj + i), given tj and tp' n (tj)). 

• evaluate ip^ (tj), . . . , ipn" 1 ^ (tj) y i a the recurrence relation (|4"T)) in Section |3~T1 (in double preci- 
sion calculations, M = 30; in extended precision calculations, M — 60). 



28 



use 20 steps of Runge-Kutta method (see Section lS.4.3[) , to solve the ODE (|5Tj) with the initial 
condition 

in the interval [w ■ (j — 1/2), 7r- (j + 1/2)] (see (|82p). The rightmost value tj+i of the solution 
is a low-order approximation of tj+i. 

• compute tj+i via Newton's method (see Section l3.4.ip . using tj+i as the initial approxima- 
tion to tj+x- On each Newton iteration, we evaluate if) n and ip' n via Taylor's method (see 
Section f3. 4. 21) . The Taylor expansion of appropriate order M about tj is used, i.e. 

Mt) = E • ( f - + ((* - ■ (we) 

fe=0 

• evaluate ip' n (tj + i) via Taylor's method. The Taylor expansion of order M — 1 is used, i.e. 

M-l ,(fc+l)/. \ 

= E )□ — • - + O (fe+i - ^) M ) • (147) 

fe=0 

In both (|146p and (|147[) . we set M — 30 for double precision calculations, and M = 60 for 
extended precision calculations. 

Observation 3. The point tj+\ approximates tj + i to at least three decimal digits (see Section l3.4.3l) . 
Subsequently, only several Newton iterations are required to obtain tj+i to essentially machine 
precision (see Observation 1 above, and also [5] for more details). Thus the cost of Step 3 is O(l) 
operations, for every integer n/2 < j < n. 

Remark 26. Obviously, on each Newton iteration one can evaluate ip n and ip' n via (|133p . (I134|) 
in Section 15. Jl rather than via (|146p . (|147p . However, this would increase the cost of each such 
evaluation from O(l) to 0(n), and the total cost of the procedure from 0{n) to 0(n 2 ) (see Remark \20\ 
in Section\5.1\). 



Step 4 (evaluation of tj and ip' n (tj) for all j < n/2). Step 3 is repeated for every integer 
n/2 < j < n. To evaluate tj and ip' n (tj) for — 1 < tj < 0, we use the symmetry of ip n about zero 
(see Theorem[T]in Section l3~Tj) . More specifically, for every integer 1 < j < n/2, we compute tj and 
ib' n {tj), respectively, via the formulae 

tj = tn+l-j (148) 

and 

^fe) = (-l)" +1 (149) 

Summary (evaluation of tj and ip' n (tj), for all j = 1, . . . , n). To summarize, the procedure 
for the evaluation of all roots of ip n in (—1,1) (as well as the derivative of ip n at these roots) is as 
follows: 

• Evaluate i m i n defined via (|143p . (|144[) (see Step 1). Cost: 0(n) operations. 

• Evaluate ip' n (t m i n ) (see Step 2). Cost: 0(n) operations. 
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• For every integer n/2 < j < n, evaluate tj+i and ip' n (tj+i) (see Step 3). Cost: 0(n) operations. 

• For every integer 1 < j < n/2, evaluate tj and tp' n (tj) (see Step 4). Cost: 0(n) operations. 

Remark 27. We observe that the algorithm described in this section not only computes the roots 
t\,...,t n of ip n in (—1,1), but also evaluates ip' n at all these roots. The total cost of this algo- 
rithm is 0(n) operations, and all the quantities are evaluated essentially to machine precision (see 
Observations 1,2,3 above). 

(n) (n) 

Remark 28. The algorithm described in this section uses the quantities \n and /3q , (i\ , . . . com- 
puted via the procedure of Section \5. 11 If n < 2c/n, then these quantities are obtained at the cost 
of 0(n + clog(c)) operations; if n > 2c/ it, then these quantities are obtained at the cost of 0(n) 
operations (see Remarks ] 161 \2(A in Section \5.1)) . 

Remark 29. As a by-product of the algorithm described in this section, we obtain a table of all 
the derivatives of tp n up to order M at all roots of tj) n in (—1,1) (here M = 30 in double precision 
calculation, and M = 60 in extended precision calculations). In other words, ip^(tj) are calculated 
for every k = 1, . . . , M and every j = 1, . . . ,n (see Step 3 above). This table can be used to evaluate 
il) n {x),il)' n {x) at an arbitrary point t\ < x < t n to essentially machine precision in 0(1) operations 
via interpolation, using the formulae (|146l) . (|147l) (see also Remark \21\ in Section \5.1\) . 



5.4 Evaluation of the Quadrature Weights 

Suppose now that n > is an integer, and that the quadrature rule S n is defined via (jlOll) in 
Section @] In this subsection, we describe an algorithm for the evaluation of the weights Wx, . . . , W n 
of this quadrature rule (see (|100j) in Section @| . The results of this subsection are illustrated in 
Table [5] and in Figure [5] (see Experiment 4 in Section 17. 2[) . 

In the description of the algorithms below, we assume that the coefficients /3q , (3^ , . . . (defined 
via (|56p in Section I3.2j) have already been evaluated (for example, by the algorithm in Section 15.11) . 
In addition, we assume that the quadrature nodes t%, . . . ,t n as well as ^' n {ti), . . . , ip' n (t n ) have also 
been computed (for example, by the algorithm of Section [5.3|) . 

An obvious way to compute W\ , . . . , W n is to evaluate (|100[) numerically. However, due to , 
the integrand (fj in (llOOp has n — 1 roots in (—1, 1), for every j = 1, . . . ,n. In particular, such 
approach is unlikely to require less that 0(n 2 ) operations. 

Rather than computing (|100[) directly, we evaluate W\ , W n by using the results of Section POl 
In the rest of this subsection, we describe two such algorithms; both evaluate W%, . . . , W n essentially 
to machine precision. One of these algorithms (based on Theorem [17]) is fairly straightforward; 
however, its cost is 0{n 2 ) operations. The other algorithm (based on Theorem [T8]) . while still rather 
simple, is also computationally efficient: its cost is 0{n) operations. 

Algorithm 1: evaluation of W±, . . . , W n in 0(n 2 ) operations. Suppose that the integer N is 
defined via (|127[) in Section 15.11 For every integer j = 1, . . . ,n, we compute an approximation Wj 
to Wj via the formula 

_ 9 2W 2N 

Wj = --_5>< n) • Q k (t,) = -—- Y.^ ■ ' VkTT/2, (150) 

Vn\ j) k=Q Vn\ j) fc=Q 

where Qk(t) and are defined, respectively, via (|68|). ((69)) and (IBTl) in Section [3~2l We observe 
that (|150[) is obtained from the identity (111 8[) in Theorem [T71 in Section |4~31 by truncating the infinite 
series at 2N terms. 
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Remark 30. Due to the combination of Remarks [ffi [77J1 in Section \3.2\ Remark\T^in Section \5.1[ 
(|127p and Theorem \17\ each weight Wj is evaluated via (|150|) essentially to machine precision (see 
also Experiment 4 in Section \7.2fy . 

Remark 31. Due to the combination of Remark \10\ in Section \3.2\ and (|127[) in Section \5.1l the 
overall cost of computing W\, . . . , W n via (|150[) is 0(n 2 ) operations. 



Algorithm 2: evaluation of Wi, . . . ,W n in 0(n) operations. This algorithm is somewhat 
similar to the procedure for the evaluation of the roots of ip n in (—1, 1) described in Section 1531 
Suppose first that t m i n is the minimal root of tp n in [0, 1). In other words, 

. _ J *(n+i)/2 = if n is odd, 

train — \ . 

{t(n+2)/2 > if n is even 

(see (I143p . (|144[) in Section [53)) . Suppose also that the function $„ : (—1,1) — > K is defined via 
(TTTTj) in Theorem ITH in Section [ 



Step 1 (evaluation of $ n (£ m i n ) and $' n (i m in))- We evaluate $ n (i m in) and 3>^(imin) via the 
formulae 

2N 2N 

$n(imin) = E 4^ ' Q*(*mln) = E ' Qfe(*min) ■ ^A: + 1/2 (152) 

fc=0 k=0 

and 

= E 4 n) • OfcCw) = E 4 n) • QUUin) • Vk+T/2, (153) 

fc=0 fe=0 

respectively (see fl!50[) in the description of Algorithm 1 above). Observe that (|152[) . f|153[) are 
obtained from the infinite expansion (|117[) in Theorem 1171 by truncation. 

Remark 32. Due to Remarks \3(A \31l the cost of Step 1 is 0(n) operations; moreover, $ n (t m i n ) 
and $^ (train) are evaluated via (|152[) . (|153|) essentially to machine precision. 

We evaluate $„ at all but the last four remaining roots of ?/>„ in [0, 1) as follows. Suppose that 
n/2 < j < n is an integer, and both 5> n (tj) and have already been evaluated. 



Step 2 (evaluation of <E> n (t, + i) and <&' n (tj + x), given & n (tj) and $' n (tj)). 

————— ————— — f2) ~ 

• use the recurrence relation ([T2D]) . (TT2Tj) (see TheoremHSlin Scction[0]) to evaluate $Jr fe), • ■ • , $ 
(here M = 60 in double precision calculations, and M — 120 in extended precision calcula- 
tions). 

• evaluate <& n (tj +1 ) via Taylor's method (see Section [233). The Taylor expansion of appropriate 
order M is used, i.e. 

M 3> {k) (t ) 

= E ^^r 2 • fe+i - ^ + (fe+i - ^) M+1 ) ( 154 ) 

(compare to (1146)) in Section I5"3j) . 
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• evaluate <b' n (tj + i) via Taylor's method. The Taylor expansion of order M — 1 is used, i.e. 

m-i $0+1)/, \ 

<(t j+i ) = - kl • - ^) fc + - *i) M ) ( 155 ) 

(compare to (11471) in Section [53]) . In both f|154[) and (|155[) , we set M = 60 for double precision 
calculations and M — 120 for extended precision calculations. 

Remark 33. For each j, the cost of Step 2 is 0(1) operations (i.e. does not depend on n). Also, it 
turns out that <$> n (tj) and Q' n (tj) are evaluated via (|154p . (|155[) respectively, essentially to machine 
precision (compare to (|146[) . ([1471) in Section \ 5.3\) . For a detailed discussion of the accuracy and 
stability of this step, the reader is referred to fSjj. 

Step 3 (evaluation of 3>„(£j) for n — 3 < j < n). For j = n — 3, n — 2, n — 1, n, we evaluate 
3>n(ij) via the formula 

= ^ o> } • Q fe (£,) = 4 n) • Q fe fe) • VfcTl72 (156) 

fc=0 fc=0 

(as in (|152[) in Step 1; see also (|150[) in the description of Algorithm 1 above). 

Remark 34. We compute $ n af the last four nodes via (|156[) rather than (1154[) . since i/ie accuracy 
of the latter deteriorates when tj is too close to 1 (interestingly, the evaluation ofip n (tj) via p46[) in 
Section \5.S\ for any j — l,...,n does not have this unpleasant feature) . Since this approach works in 
practice, is cheap in terms of the number of operations and eliminates the accuracy problem, there 
was no need in a detailed analysis of the issue (see, however, J^j for more details). 

Step 4 (evaluation of & n (tj) for 1 < j < n/2). Due to the combination of Theorem [T71 in 
Section 14.31 and (|6"5|) in Section 13.21 the function is symmetric about the origin. We use this 
observation to evaluate <&„ (tj ) via the formula 

- (-1)" +1 ■ $„((„ +H ), (157) 

for every j — 1,2, ... , n/2. 



Step 5 (evaluation of Wi, . . . , W n ). For every j 
to Wj from & n (tj) and ip' n (tj) via the formula 



1, . . . , n, we compute an approximation Wj 



(see (fHgJl in Theorem Q7| in Section O . 



Remark 35. Due to the combination of Remarks HIR 
essentially to machine precision. This algorithms requires 0(n) operations (compare to Remark \31)) . 



(158) 



Algorithm 2 evaluates all Wi, . ■ ■ , W n 



Remark 36. Algorithm 2 described in this section uses some of the quantities evaluated by the proce- 
dures of Sections \5.1[ \5.S\ If n < 2c/n, then the cost of obtaining these quantities is O (n + clog(c)) 
operations; if n > 2c/tt, then the cost of obtaining these quantities is 0(n) operations (see Re- 
marks\F^\M in Section^). 
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6 Numerical Results 



In this section, we demonstrate the performance of the quadrature rules from Section [U All the 
calculations were implemented in FORTRAN (the Lahey 95 LINUX version), and carried out in 
double precision. Extended precision calculations were used for comparison and verification (in 
extended precision, the floating point numbers are 128 bits long, as opposed to 64 bits in double 
precision). 



5 n (e ) in double precision 
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Figure 1: The quadrature error vs |A n |, with c — 1000 and n — 682. _ffere A„ = -.60352E-15. 
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Experiment 1. Here we demonstrate the performance of the quadrature rule S n (see (jlOip in 
Section @| on exponential functions. We proceed as follows. We choose, more or less arbitrarily, 
the band limit c and the prolate index n. Next, we evaluate the quadrature nodes t\,...,t n and 
the quadrature weights W%, . . . ,W n via the algorithms of Sections 15.31 15.41 respectively. Also, we 
evaluate |A„| via the algorithm in Section T5.2I Then, we choose a real number < a < 2, and 
evaluate the integral of e %cax over — 1 < x < 1 via the formula 

C e iacx dx = f 1 cos(acx) dx = 2sm ( ac \ ( 159 ) 
J-i J-i ac 

Also, we use S n to approximate (| 159[) via the formula 

/l n 
e lacx dx^Y^ e lcat > ■ Wj (160) 

1 3=1 
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5 n (e ) in extended precision 




0.5 1 1.5 



a 

Figure 2: The quadrature error vs |A„|, with c = 1000 and n = 682. Here X n — -.60352E-15. 



A (n), A (n) vs. |A, |, in double precision 




Figure 3: The maximal quadrature errors Ai(n), Aa(n) vs \X n \, with c — 1000. 




Figure 4: The maximal quadrature errors Ai(n), Aa(n) us |A„|, irai/i c = 1000. 



(see (|102p in Section 0]). Finally, we evaluate the quadrature error 8 n (e lacx ) via the formula 

2 sin(ac) 



5 n (e iacx ) 



V" e lcat > ■ Wj 



3=1 



(161) 



(see (fTU5|) in Section g]). 

In Figure [TJ we display the results of this experiment. The band limit and the prolate index were 
chosen to be, respectively, c = 1000 and n = 682. For this choice of parameters, A„ = -.60352E-15. 
In this figure, we plot the quadrature error (|161[) as a function of the real parameter a, for < a < 2, 
on the logarithmic scale. The calculations are carried out in double precision. 

We make the following observations from Figure [TJ The quadrature error is essentially zero up 
to machine precision e, for all real < a < 2. In other words, for this choice of parameters, the 
quadrature rule S n integrates the functions of the form f(x) = e %cax with < a < 1 exactly, for all 
practical purposes. It is perhaps surprising, however, that such functions are integrated exactly via 
S n even when 1 < a < 2. In other words, the quadrature rule S n (corresponding to band limit c 
and |A„| rs e) integrates exactly the exponential functions with the band limit up to 2c. 

To get a clearer picture, we repeat this experiment in extended precision. In Figure (2[ we plot 
the quadrature error (|161[) as a function of the real parameter a, for < a < 2, on the logarithmic 
scale. In other words, Figure [5] is a version of Figure [1] in extended precision. 

We make the following observations from Figure If < a < 1, then the quadrature rule S n 
integrates the functions of the form f(x) — e tcax up to the error of order |A n | 2 (in Figured] we used 
double precision calculations and thus did not have enough digits to see this phenomenon). On the 
other hand, for 1 < a < 2 the quadrature rule S n integrates such functions up to the error roughly 
|A„|. In other words, the quadrature rule S n (corresponding to band limit c and |A„| we) integrates 
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the functions of band limit up to c up to e 2 (rather than e); on the other hand, the functions of 
band limit between c and 2c are integrated up to e. 

Explanation. These observations admit the following (somewhat imprecise) explanation (see 
[2"5] . [2"S] for more details). Suppose that a > is a real number. Due to (|27D and Theorem Q] in 
Section 13.11 



eiacx = J- \ m il> m {a)i> m {x), (162) 

m=0 

for all real — 1 < x < 1. Moreover, 

f 1 e »- da; = 2sm(ac) = g ^ (o) ^ (0)< (163) 

We combine ([TCTl) . (fl6"2l . (|T6^1) to obtain 

r, • / \ n oo / n \ 

- W i ■ e " Qtj = E A ™^™(«) A m ^ m (0) - J2 W^ m {tj) ■ (164) 

Obviously, the quadrature error 5 n {ipm) (see (|173p ) is zero for odd m. Also, 5 n (ip m ) rapidly increases 
as a function of even < m < n; moreover, 6 n (?p m ) is of order |A„| when m < n is an even integer 
close to n (see Conjectures [H [3] in Section I7TT1 and Theorem IT41 in Section |4TT|) . Therefore, roughly 
speaking, 

n— 1 / n \ 

E A m ^m(a) A ro Vm(0) - 51 « l*n| a ■ ^-i(o). (165) 

m=0 \ j = l / 

On the other hand, due to the fast decay of |A m | (see Theorems ® [7] in Section ROT) . 

00/ n \ 

^ A m ^m(a) A m Vm(Q) - E « |A„| 2 . (166) 

m=n y j=i y 

Finally, the following approximate formula appears in |25j . |26) . in a slightly different form: suppose 
that n > is an integer, that Xn > c 2 , and that < a < 2 is a real number. Then, 

\Ma)\ = [° { f]\, (167) 
IW ;l [OdA^- 1 ), Ka<2. V ^ 

It follows from the combination of (|165[) , (|166[) , (|167[) that the quadrature error (|161[) is expected to 
be of the order |A„| 2 • ^/n, if < a < 1. On the other hand, the quadrature error (I161[) is expected 
to be of the order |A n |, if 1 < a < 2. Figures [TJ [3J [3J IH support these somewhat vague conclusions. 

We summarize this crude analysis, supported by the observations above, in the following conjec- 
ture about the quadrature error (|161|) for < a < 2. 

Conjecture 1. Suppose that c > and a > are real numbers, and that n > 2c/7r is an integer. 
Suppose also that 5 n (e lcax ) is defined via (| 103[) in Definition^ in Section^ IfO<a<l, then 

/i n 
e lcax dx - V e lcat > ■ Wj « |A„| 2 • V7{, (168) 
1 i=i 
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where X n is that of (|27[) in Section \3.1l If on the other hand, 1 < a < 2, then 



S n (e* cax ) 



/l n 
e lcax dx-Y] 

1 3=1 



e lcat > ■ Wi 



|A n |. 



(169) 



We repeat the above experiment with various values of n, and plot the results in Figure [3J This 
figure also corresponds to band limit c = 1000. We plot the following three quantities as functions 
of the prolate index n that varies between 637 ~ 2c/ir and 700. First, we plot |A n |. Second, we plot 
the maximal quadrature error Ai(n) defined via the formula 



Ai(n) 



max 5 n (e l 

0<a<l 



max 

0<a<l 



2 sin(ac) 



w. 



(n) 



(170) 



where t^\ . . . , t„ 1 ' and W^'* 1 , ■ ■ . , Wn""' are, respectively, the notes and weights of the quadrature 
rule S n (see (|10ip in Section [4]). Finally, we plot the maximal quadrature error A 2 (n) defined via 
the formula 



r(n) 



A 2 (n) 



max 8 n (e l 

Ka<2 



max 

Ko<2 



2 sin(ac) 



Eicat 



00 . w {n) 

3 



(171) 



We observe that in (|170p the parameter a varies between and 1, and in (I171[) the parameter a 
varies between 1 and 2. In other words, Ai(n) is the maximal quadrature errors of S n for the 
exponential functions of band limits up to c, and A2(n) is the maximal quadrature error of S n for 
the exponential functions of band limit between c and 2c. 

We make the following observations from Figure[3] As long as |A„| is less than roughly 10~ 7 « \fe 
(with e the machine precision), Ai(n) is roughly equal to |A„| 2 . On the other hand, Ai(n) is zero 
up to machine precision once |A n | > 10 -7 . These observations are in agreement with Conjecture Q] 
above. 

We also observe that A 2 (n) is roughly of order |A„|, as long as |A„| > e. On the other hand, 
when A„ is zero to machine precision, so is A2(n) (see Conjecture [T]). 

We repeat this experiment in extended precision, and plot the results in Figure 2] In other 
words, Figure [4] is a version of Figure [3] in extended precision. We observe the same phenomenon: 
Ai(n) is of order |A n | 2 , and A 2 (n) is of order |A„| (as long as we do not run out of digits to see it; 
if, for example, |A„| is below the machine zero so are both Ai(n) and ^{n)). In other words, the 
quadrature error of S n for exponential functions with band limit up to c is of order |A ra | 2 , and the 
quadrature error of S n for exponential functions with band limit between c and 2c is of order |A n |, 
which supports Conjecture [T] 



7 Numerical Illustration of Analysis in Section [4] 

In this section, we illustrate the analytical results from Section 2] and the performance of the algo- 
rithms described in Section [5] All the calculations were implemented in FORTRAN (the Lahey 95 
LINUX version) , and carried out in double precision. Extended precision calculations were used for 
comparison and verification (in extended precision, the floating point numbers are 128 bits long, as 
opposed to 64 bits in double precision). 



37 



7.1 Quadrature Error and its Relation to |A n | 

In this section, we describe several numerical experiments that illustrate the quadrature error (see 
(TTUTjl . (TiTj3")l in Section El) and its relation to |A„|. 



m 


AmV'm(O) 


S n (ipm), double precision 


Sniipm), extended precision 





0.70669E+00 


0.44409E-15 


0.33258E-26 


2 


0.49581E+00 


0.16653E-15 


0.22426E-25 


4 


0.42581E+00 


0.13323E-14 


0.26756E-23 


6 


0.38527E+00 


0.21649E-14 


0.19692E-21 


8 


0.35695E+00 


0.22760E-14 


0.91546E-20 


10 


0.33516E+00 


0.16653E-14 


0.29148E-18 


12 


0.31730E+00 


0.23870E-14 


0.88165E-17 


14 


0.30201E+00 


0.24980E-14 


0.21007E-15 


16 


0.28844E+00 


0.11102E-14 


0.35574E-14 


18 


0.27604E+00 


0.59230E-13 


0.57028E-13 


20 


0.26435E+00 


0.83716E-12 


0.83954E-12 


22 


0.25299E+00 


0.89038E-11 


0.89011E-11 


24 


0.24150E+00 


0.76862E-10 


0.76864E-10 


26 


0.22919E+00 


0.65870E-09 


0.65870E-09 


28 


0.21377E+00 


0.45239E-08 


0.45239E-08 


30 


0.18075E+00 


0.19826E-07 


0.19826E-07 


32 


0.10038E+00 


0.68548E-07 


0.68548E-07 


34 


0.27988E-01 


0.33810E-06 


0.33810E-06 


36 


0.49822E-02 


0.27232E-05 


0.27232E-05 


38 


0.70503E-03 


0.22754E-04 


0.22754E-04 



Table 1: Illustration of Theorem \14\ with 
0.12915E-03. See Experiment 2. 



50 and n = 40. For these parameters, A r 



Experiment 2. Here we illustrate Theorem [T4]in Section l4~T1 We choose, more or less arbitrarily, 
band limit c and prolate index n. We evaluate A„ and the quadrature rule S n defined via (|10ip 
in Section @] via the algorithms of Sections 15.11 15.21 15.31 15.41 respectively. Then, we choose an even 
integer < m < n, and evaluate A m , V'm(O), and ip m {tj) for all j = 1, . . . , n, via the algorithms of 
Sections 15.11 15.21 All the calculations are carried out in double precision. 

We display the results of this experiment in Table Q] The data in this table correspond to c = 50 
and n = 40. Table [1] has the following structure. The first column contains the even integer m, that 
varies between and n — 2. The second column contains X m Tp m (0) (we observe that 



A m ^m(0) 



^m(t) dt, 



due to (|2~T)) in Section [OJ- The third column contains the quadrature error 



(see f|103[) in Section 0]), computed in double precision. 



A m -0m(O) - 2J V>m(*i) • Wj 
3=1 



(172) 



(173) 



38 



Then, we repeat all the calculations in extended precision; the last column of Table [T] contains 
Sntym) defined via (|173[) (the same quantity as in the third column evaluated in extended precision). 

We make the following observations from Table [TJ We note that A m , m (O) is always positive 
and monotonically decreases with m. We also note that 5 n (ip m ) (evaluated in double precision) is 
close to the machine accuracy for small m, and grows up to w 2 ■ 10~ 5 for m = 38. Also, 5 n (ifj m ) 
is bounded by |A„|, for all values of m in Table [T] (in this case, |A„| = 0.12915E-03). Finally, 
^n(V'm) (evaluated in extended precision) is a monotonically increasing function of even < m < n 
(obviously, 5 n (ip m ) — for odd m). 

We summarize these observations in the following conjecture. We have not fully investigated the 
phenomenon described in this conjecture; see, however, Theorem fT?] in Section I4.ll Conjecture |3] 
below, Figure [5] and Table [3] (see also [25], [26] for additional details and analysis). 

Conjecture 2. Suppose that c > 1 is a real number, that n > 2c/n is an integer, and that the 
quadrature rule S n is defined via (|101[) in Section^ Then, the quadrature error 5 n (ip m ) defined 
via (I173P above is a monotonically increasing function of even < m < n. Moreover, in double 
precision calculations 5 n (ip m ) is zero up to machine precision for all < m < 2c/ir. 



In (|106p in Theorem 1141 we provide an upper bound on 5 n (ip r , 
on m; more specifically, for every m = 0, . . . ,n — 1, 



This bound does not depend 



0~n{4>m) 



A m ^> m (0) - y^tpm(tj 



3 = 1 



■■ IA„I- I 24 • log | — ) ~0-\, 



(174) 



On the other hand, according to Table[T]the quadrature error 5 n (?p m ) is bounded by |A n | alone, for 
all even < m < n (obviously, Snitpm) = for all odd m). 

In Figure [S] we display the results of the same experiment with different choice of parameters c 
and n. Namely, we choose c = 10000 and plot A m V'm(0) as a function of even < m < 6425, on the 
logarithmic scale (solid line). In addition, we plot the quadrature error S n (ip m ) as a function of m, 
for four different values of n: n — 6393 (dashed line), n — 6401 (circles), n = 6414 (triangles), and 
n = 6425 (pluses). The corresponding values of |A n | are displayed in Table [2] 



n 


6393 


6401 


6414 


6425 


|An| 


0.43299E-07 


0.54119E-09 


0.33602E-12 


0.52616E-15 



Table 2: Values of \X n \ for c — 10000 and different choices of n. 

We make the following observations from Figure [S] First, the quantities A rra '(/' m (0) are of the 
same order of magnitude for all m < 2c/ir, and decay rapidly with m for m > 2c/tt. Also, for each 
value of n, the quadrature error 8 n (ip m ) is essentially zero for all m < 2c/7r and increases rapidly 
with m for m > 2c/n. Nevertheless, 5 n (jp m ) is always bounded from above by |A„|, for each n and 
each m < n. See also Tables Q] |3] and Conjecture |3] below. 

We repeat the experiment above with several other values of band limit c and prolate index n. 
The results are displayed in Table [3] This table has the following structure. The first and second 
column contain, respectively, the band limit c and the prolate index n. The third column contains 
the even integer < m < n (the values of m were chosen to be close to n). The fourth column 
contains A m '0 m (O). The fifth column contains (|173[) . The last column contains |A„|. 

We make the following observations from Table [3] First, for each of the seven values of c, the 
three indices n were chosen in such a way that |A n | is between 10 -12 and 10~ 7 . The values of the 
band limit c vary between 250 (the first three rows) and 16000 (the last three rows). For each 
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Quadrature error 



10" 



10" 



10" 



10" 




Figure 5: The quadrature error <5 n (Vv> 



/l n 
^j m {t) dt-^^ m {t )-W 3 

- 1 3 = 1 



as a function of even 



m < n, for four different values of n and c = 10000, vs. A m "i/> m (0). See Experiment 2. 



40 



c 


n 


m 


A m Vm(0) 


f 1 -A 

J-l j=l 


|An| 


250 


179 


178 


0.28699E-07 


-.52496E-08 


0.18854E-07 


250 


184 


182 


0.68573E-09 


-.38341E-10 


0.16130E-09 


250 


188 


186 


0.14108E-10 


-.68758E-12 


0.30500E-11 


500 


339 


338 


0.52368E-07 


-.13473E-07 


0.40938E-07 


500 


345 


344 


0.37412E-09 


-.86136E-10 


0.27418E-09 


500 


350 


348 


0.12148E-10 


-.99816E-12 


0.35537E-11 


1000 


659 


658 


0.42709E-07 


-.14354E-07 


0.38241E-07 


1000 


665 


664 


0.51665E-09 


-.15924E-09 


0.43991E-09 


1000 


671 


670 


0.52494E-11 


-.15024E-11 


0.42815E-11 


2000 


1297 


1296 


0.41418E-07 


-.17547E-07 


0.41740E-07 


2000 


1304 


1302 


0.77185E-09 


-.15036E-09 


0.37721E-09 


2000 


1311 


1310 


0.31078E-11 


-.11386E-11 


0.28754E-11 


4000 


2572 


2570 


0.54840E-07 


-.15493E-07 


0.33682E-07 


4000 


2579 


2578 


0.43032E-09 


-.20771E-09 


0.46141E-09 


4000 


2587 


2586 


0.28193E-11 


-.12805E-11 


0.29164E-11 


8000 


5119 


5118 


0.43268E-07 


-.26751E-07 


0.52899E-07 


8000 


5128 


5126 


0.50230E-09 


-.16395E-09 


0.33442E-09 


8000 


5136 


5134 


0.50508E-11 


-.15448E-11 


0.32132E-11 


16000 


10213 


10212 


0.42725E-07 


-.30880E-07 


0.56568E-07 


16000 


10222 


10220 


0.69663E-09 


-.28201E-09 


0.52821E-09 


16000 


10231 


10230 


0.34472E-11 


-.22162E-11 


0.42902E-11 



Tabic 3: Relation between the quadrature error and \X n \. See Experiment 2. 
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n, the value of m is chosen to be the largest even integer below n. This choice of m yields the 
smallest X m ip m (0) and the largest quadrature error 8 n {tpm) among all m < n (see also Table [JJ and 
Figure [SJ . Obviously, for this choice of m the eigenvalues X m and A„ are roughly of the same order 
of magnitude. Wc also observe that for all the values of c,n,m, the quadrature error 5 n (ip m ) is 
bounded from above by |A„| (and is roughly equal to |A„|/2). In other words, the upper bound on 
Sn(fpm) provided by Theorem [Jj] (see (|174p ) is somewhat overcautious. 
We summarize these observations in the following conjecture. 

Conjecture 3. Suppose that c > is a positive real number, and that n > 2c/7r is an integer. 
Suppose also that < m < n is an integer. Suppose furthermore that (^(VVn) is defined via (|103|) in 
Definition [H in Section [7J Then, 



/l n 
1p m {s) ds-y^ i>m{tj) ■ W. 

- 1 S=l 



< |A„ 



(175) 



where X n is that of (|27p in Section \ c 3.1\ 



Remark 37. The ineguality (|175p in Conjecture [3] is stronger than the ineguality (|106[) in Theo- 
rem \14\ On the other hand, as opposed to Theorem \14\ Conjecture [21 has been only supported by 
numerical evidence. 



Experiment 3. Here we illustrate Theorems ITSl [TBI in Section l4~2l We proceed as follows. We 
choose, more or less arbitrarily, the band limit c > and the accuracy parameter e > 0. Then, we 
use the algorithm of Section 15.21 to find the minimal integer m such that |A m | < e. In other words, 
we define the integer n\{e) via the formula 

m(e) = min{m > : |A m |<e}. (176) 

Also, we find the minimal integer such that the right-hand side of (|106[) in Theorem [TJ] in Section I4TT1 
is less that e. In other words, we define the integer 712(e) via the formula 

na(e) = min jm > : |A m | • (^24 • log f p^J + 6 • Xmj < £ j • (177) 

Next, we evaluate the integer n^(e) via the formula (|110[) in Theorem 1151 In other words, 

/2c a(e) , / 16ec \\ . m . 

Me) = floor (- + ^- log ( wy )) (178) 

where a(e) is defined via (|109[) in Theorem [T5l Finally, we evaluate the integer n^e) via the 
right-hand side of (II 15|) in Theorem [TBI In other words, 

n 4 ( £ )=floor('^+('lO+^log(c) + i-logi) - log (5)) ■ (179) 

In both (|178p and (|179p . floor(a) denotes the integer part of a real number a. 

We display the results of this experiment in Table @] This table has the following structure. 
The first column contains the band limit c. The second column contains the accuracy parameter e. 
The third column contains n\(e) defined via (|176p . The fourth column contains ^(e) defined via 
(|177p . The fifth column contains 713(e) defined via (|178p . The sixth column contains 714(e) defined 
via (|179l) . The seventh column contains |A ni ( e )|. The last column contains |A„ 2 ( e )|. 
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c 


e 


m(e) 


n 2 (e) 


n 3 (e) 


n 4 (e) 


\Ki(e)\ 


\K 2 (e)\ 


250 


in- 
1U 


10 


184 


198 


277 


303 


0.60576E-10 


0.86791E-16 


250 




25 


216 


227 


326 


386 


0.31798E-25 


0.14863E-30 


250 


1U 


50 


260 


270 


393 


525 


0.28910E-50 


0.75155E-56 


500 


in- 
1U 


10 


346 


362 


460 


488 


0.49076E-10 


0.60092E-16 


500 


1U 


25 


382 


397 


520 


583 


0.54529E-25 


0.19622E-31 


500 


1 n~ 


50 


433 


446 


607 


742 


0.82391E-50 


0.38217E-56 


1000 


in- 
1U 


10 


666 


687 


803 


834 


0.95582E-10 


0.92947E-17 


1000 


1 0" 


25 


707 


725 


875 


942 


0.97844E-25 


0.14241E-31 


1000 


1U 


50 


767 


783 


981 


1120 


0.39772E-50 


0.56698E-57 


2000 


1U 


10 


1305 


1330 


1467 


1500 


0.95177E-10 


0.25349E-17 


2000 


10" 


25 


1351 


1373 


1550 


1619 


0.86694E-25 


0.27321E-32 


2000 


io- 


50 


1418 


1438 


1675 


1818 


0.88841E-50 


0.22795E-57 


4000 


10" 


10 


2581 


2610 


2768 


2804 


0.70386E-10 


0.64396E-18 


4000 


io- 


25 


2632 


2658 


2862 


2935 


0.57213E-25 


0.53827E-33 


4000 


io- 


50 


2707 


2730 


3007 


3154 


0.56712E-50 


0.88819E-58 


8000 


io- 


10 


5130 


5163 


5344 


5383 


0.59447E-10 


0.22821E-18 


8000 


io- 


25 


5185 


5216 


5450 


5526 


0.87242E-25 


0.16237E-33 


8000 


io- 


50 


5268 


5296 


5614 


5765 


0.95784E-50 


0.23927E-58 


16000 


io- 


10 


10225 


10264 


10468 


10509 


0.63183E-10 


0.37516E-19 


16000 


io- 


25 


10285 


10321 


10585 


10664 


0.85910E-25 


0.41416E-34 


16000 


io- 


50 


10377 


10409 


10769 


10923 


0.51912E-50 


0.56250E-59 


32000 


io- 


10 


20413 


20457 


20686 


20730 


0.62113E-10 


0.12818E-19 


32000 


io- 


25 


20478 


20519 


20815 


20897 


0.78699E-25 


0.12197E-34 


32000 


io- 


50 


20577 


20615 


21018 


21176 


0.96802E-50 


0.15816E-59 


64000 


io- 


10 


40786 


40837 


41092 


41139 


0.89344E-10 


0.28169E-20 


64000 


io- 


25 


40857 


40903 


41232 


41318 


0.66605E-25 


0.39212E-35 


64000 


io- 


50 


40964 


41008 


41454 


41616 


0.85451E-50 


0.28036E-60 


10 s 


io- 


10 


636669 


636747 


637115 


637174 


0.79326E-10 


0.13385E-22 


10 6 


io- 


25 


636759 


636832 


637301 


637400 


0.77413E-25 


0.15758E-37 


10 6 


io- 


50 


636899 


636968 


637600 


637778 


0.69235E-50 


0.15801E-62 



Table 4: Illustration of Theorems ] 151 [751 See Experiment 3. 
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Suppose that c > is a band limit, and n > is an integer. We define the real number Q(c,n) 
via the formula 



Q(c, n) = max I 6 n (ip m ) 



/l n 
1> m {t)dt-^i/i m {tj)-W. 

- 1 3 = 1 



0<m<n-l \ , (180) 



where t%, . . . , t n and W\, . . . , W n are defined, respectively, via (l98l) . (jlOOl) in Definition [2] in Section[4] 
In other words, Q{c, n) is the maximal error to which the quadrature rule S n defined via (jlOip 
integrates the first n PSWFs. 

We make the following observations from TableHJ We observe that Q(c,ni(e)) < e, due to the 
combination of Conjecture [3] in Section 17711 and (|176D . (|180[) . In other words, numerical evidence 
suggests that the quadrature rule 5 ni ( e ) integrates the first ni(e) PSWFs up to an error less than 
e (see Remark [37]). On the other hand, we combine Theorem [Ml in Section |4~T1 with (|177p . fll80[) . 
to conclude that the quadrature rule S n2 i e \ has been rigorously proven to integrate the first 772(e) 
PSWFs up to an error less than e. In both Theorem [T~41 and Conjecture^ we establish upper bounds 
on Q(c, n) in terms of |A„|. The ratio of |A ni ( £ )| to |A„ 2 ( e )| is quite large: from about 10 6 for c = 250 
and e = 10" 10 , 10~ 25 , 10" 50 (see the first three rows in Table SJ), to about 10 10 for c = 64000 and 
e = lO" 10 , 10" 25 , 10~ 50 , to about 5 • 10 12 for c = 10 6 and e = 10~ 10 , 10" 25 , 10~ 50 , (see the last six 
rows in Table 2]). On the other hand, the difference between 772(e) and fii(e) is fairly small; for 
example, for e = 10 -50 , this difference varies from 10 for c = 250 to 23 for c = 4000, to merely 44 
for c = 64000 and 69 for as large c as c = 10 6 . 

As opposed to n\{e) and n2(£), the integer ris(e) is computed via the explicit formula (|178[) that 
depends only on c and e (rather than on |A„| and \n, that need to be evaluated numerically); this 
formula appears in Theorem [TBI The convenience of f| 1 78[> vs. (| 1 76[) . (| 1 77[> comes at a price: for 
example, for e = 10~ 50 , the difference between n%(e) and 712(e) is equal to 123 for c = 250, to 446 for 
c = 64000, and to 632 for c = 10 6 . However, the difference 713(e) — 772(e) is rather small compared 
to c: for example, for e = 10~ 50 , this difference is roughly 4 • (log(c)) , for all values of c in TableHJ 

Furthermore, we observe that 774(e) is computed via the explicit formula (|179|) that depends only 
on c and e. This formula can be viewed as a simplified version of (|178[) (see Theorems 1151 II 6[) : in 
particular, 774(e) is greater than 773(e), for all c and e. 

We summarize these observations as follows. Suppose that the band limit c > and the accuracy 
parameter e > are given. According to Theorem 1151 for any 77 > 773(e) the quadrature rule S n 
defined via (|101[) in Section @] is guaranteed to integrate the first 77 PSWFs to precision e (see 
Definition [1] in Section l2~Tj) . On the other hand, numerical evidence (see Experiment 2) suggests 
that the choice 77 > 773(e) is overly cautious for this purpose; more specifically, S n integrates the 
first 77 PSWFs to precision e for every 77 between 771(e) and 773(e) as well. In this experiment, we 
observed that the difference between the "theoretical" bound 773(e) and "empirical" bound ni(s) is 
of order (log(c)) 2 , and, in particular, is relatively small compared to both 771(e) and 773(e) (which 
are of order c). 

Finally, we observe that 

2c 2 1 
771(e) < — + — •(logc)-log-, (181) 

for all the values of c and e in Table 2J Combined with some additional numerical experiments by 
the authors, this observation leads to the following conjecture (see also Theorem [7] in Section |3~T1 for 
a rigorously proven and more precise statement). 

Conjecture 4. Suppose that c > 1 a?7<i < e < 1 are real numbers. Suppose also that n > is an 
integer, and that 

77 > — + 10+^j-(logc)-logi. (182) 
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Then, 



|A„|<£, (183) 

where X n is that of (|27p in Section \3.1l 

7.2 Quadrature Weights 

In this section, we illustrate the results of Section 14.31 and the algorithms of Section 15.41 



3 




Wj - Wj 


f7> M / ( n +i)/2(V'„(0)) 

3 (^(* 3 )) 2 -(i-t?) 


1 


0.7602931556894E-02 


0.00000E+00 


-.55796E-11 


2 


0.1716167229714E-01 


0.00000E+00 


-.55504E-10 


3 


0.2563684665002E-01 


0.00000E+00 


-.21825E-12 


4 


0.3278512460580E-01 


0.00000E+00 


-.11959E-09 


5 


0.3863462966166E-01 


0.16653E-15 


0.82238E-11 


6 


0.4334940472363E-01 


0.22204E-15 


-.16247E-09 


7 


0.4713107235981E-01 


0.22204E-15 


0.11270E-10 


8 


0.5016785516291E-01 


0.19429E-15 


-.18720E-09 


9 


0.5261660773966E-01 


0.26368E-15 


0.10495E-10 


10 


0.5460119701692E-01 


0.29837E-15 


-.20097E-09 


11 


0.5621699326080E-01 


0.17347E-15 


0.81464E-11 


12 


0.5753664411864E-01 


0.12490E-15 


-.20866E-09 


13 


0.5861531690539E-01 


0.10408E-15 


0.55098E-11 


14 


0.5949490764741E-01 


0.23592E-15 


-.21301E-09 


15 


0.6020725336886E-01 


0.13184E-15 


0.31869E-11 


16 


0.6077650804037E-01 


0.18041E-15 


-.21545E-09 


17 


0.6122088420703E-01 


0.48572E-16 


0.14361E-11 


18 


0.6155390478472E-01 


0.83267E-16 


-.21675E-09 


19 


0.6178529976346E-01 


0.11102E-15 


0.36146E-12 


20 


0.6192162112196E-01 


0.48572E-16 


-.21732E-09 


21 


0.6196665001384E-01 


0.00000E+00 


0.00000E+00 



Table 5: Quadrature weights (|100[) with c — 40, n = 41. A rl = i0.69857E-08. See Experiment 4- 



Experiment 4. In this experiment, we choose, more or less arbitrarily, band limit c and pro- 
late index n. Then, we compute ti,...,t n (see and ip' n (ti), . . . , ip' n (t n ) via the algorithm of 
Section 15751 Also, we evaluate t/Vi(0), ip' n (0) via the algorithm of Section [S~T1 Next, compute approx- 
imations W%, . . . , W n to Wx, . . . , W n via Algorithm 1 in Section [5741 (in particular, Wj is evaluated 
via (|150[) for every j = 1, . . . , n). Also, we compute approximations W%, . . . , W n to Wi,..., W n via 
Algorithm 2 in Section 15.41 All the calculations are carried out in double precision. 

We display the results of this experiment in Table [5] The data in this table correspond to c = 40 
and n = 41. Table [5]has the following structure. The first column contains the weight index j, that 
varies between 1 and (n + l)/2 = 21. The second column contains Wj ( an approximation to Wj 
evaluated by Algorithm 2 in Section I5~4")) . The third column contains the difference between Wj and 



45 



Wj (evaluated via (|150[) by Algorithm 1). The last column contains the difference 



w< 



(n+l)/2 ■ 



(^fe)r-(i-^) 



(184) 



(see (|125l) in Remark [T4]) . 

In Figure [51 we plot the weights Wj as a function of j = 1, . . . , n. Each Wj is plotted as a circle 
above the corresponding node tj. 

We make the following observations from Table [5] First, due to the combination of Theorems IT71 
[TBI in Section l4~3l the value in the third column would be zero in exact arithmetic. We observe that, 
indeed, this value is zero up to the machine precision^which confirms Remarks 1301 1351 in Section [53] 
(We note that, for j = 1,2,3,4 and j = 21, both Wj and Wj are evaluated via (|150l) . and hence 
the value in the corresponding rows is exactly zero). In particular, either of the two approximations 
Wj , Wj can be used to evaluate Wj to essentially machine precision. 

We also observe that all of the weights Wi, . . . , W n are positive (see Theorem [T9l and Remark [TBI . 
Moreover, Wj grow monotonically as j increases to (n + l)/2. Finally, we observe that, for all 
j = 1, . . . , n, the value (I184j) in the last column is of the order |A„| (see Remark fl"4"|) . 



Quadrature weights 




Figure 6: The quadrature weights W\, . . . , W„ with c = 40, n = 41. See Experiment 4- 
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